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]';