%
% modified Euler's Method
%
% [tv,yv,n]=modeul(a,b,y0,h);
%
function [tv,yv,n]=modeul(a,b,y0,h);
n=floor((b-a)/h);
yv(1)=y0;
tv(1)=a;
for ii=2:n+1;
    f1=fun(tv(ii-1),yv(ii-1));
    tempy=yv(ii-1)+h*f1;
    f2=fun(tv(ii-1)+h,tempy);
    yv(ii)=yv(ii-1)+h*.5*(f1+f2);
    tv(ii)=tv(ii-1)+h;
end;