function fp=diff_1_2(f, x, h)   
    x_p2h = x + 2*h;
    x_ph  = x + h;
    x_mh  = x - h;
    x_m2h = x - 2*h;
    fp=(-1.0/12.0*f(x_p2h) + 2.0/3.0*f(x_ph) - 2.0/3.0*f(x_mh)+1.0/12.0*f(x_m2h))./h;
end