function [x,y]=fdiff(funcs,tspan,x0f,n)
   t0=tspan(1);tfinal=tspan(2); ga=x0f(1); gb=x0f(2);
   h=(tfinal-t0)/n;
   for i=1:n, x(i)=t0+h*(i-1); end,
   x0=x(1:n-1); t=-2+h^2*feval(funcs,x0,2); tmp=feval(funcs,x0,1);
   v=1+h*tmp/2; w=1-h*tmp/2; b=h^2*feval(funcs,x0,3);
   b(1)=b(1)-w(1)*ga; b(n-1)=b(n-1)-v(n-1)*gb; b=b'; A=diag(t);
   for i=1:n-2, A(i,i+1)=v(i); A(i+1,i)=w(i+1); end
   y=inv(A)*b;  x=[x tfinal]; y=[ga; y; gb]';