數值微分

function numerical_diff = ND(f,x,h)

f_1=f(x-2*h);
f_2=f(x-h);
f_3=f(x+h);
f_4=f(x+2*h);

numerical_diff = 1/12*(f_1-8*f_2+8*f_3-f_4)./h;

end