function y=der2p(f,x0,h)
y=1/h*(f(x0+h/2)-f(x0-h/2));