Here is the matlab/freemat code I got to solve ODE numerically using backward Euler method. However the results are inconsistent with my textbook results, and sometimes even ridiculously inconsistent. Please point out what is wrong with the code.

I have already asked this question on mathoverflow.com no help there, hope someone here can help.

function [x,y]=backEuler(f,xinit,yinit,xfinal,h)

%f - this is your y prime
%xinit - initial X
%yinit - initial Y
%xfinal - final X
%h - step size

n=(xfinal-xinit)/h; %Calculate steps

%Inititialize arrays...
%1st elements take xinit and yinit corespondigly, the rest fill with 0s
x=[xinit zeros(1,n)]; 
y=[yinit zeros(1,n)];

%Numeric routine
for i=1:n
x(i+1)=x(i)+h;
ynew=y(i)+h*(f(x(i),y(i)));
y(i+1)=y(i)+h*f(x(i+1),ynew);
end

end