clear;
f=inline('(y^2 + 2*t*y)/(3+t^2)','t','y');
yexact=inline('(3+t^2)/(6-t)','t');
hlist=[0.1 0.05 0.025 0.0125];
for hi=1:length(hlist)
    clear t y1 y2 yex;
    h=hlist(hi);
    disp(h);
    t=0.0:h:1.0;
    N=length(t);  % total number of time steps
    % initial conditions
    y1(1)=0.5;      % approximate Euler solution
    y2(1)=0.5;      % approximate Improved Euler solution
    yex(1)=0.5;     % exact solution (analytic)
    % start the approximation stepping
    for n=1:(N-1)
        y1(n+1)=y1(n) + h*f(t(n),y1(n));    % Euler's Method
        y2(n+1)=y2(n) + (h/2)*(f(t(n),y2(n)) +...
                               f(t(n+1),y2(n)+h*f(t(n),y2(n))));
        yex(n+1)=yexact(t(n+1));              % exact solution
    end
    figure;                                 % open up a new figure
    %figure(10);
    %hold on;
    plot(t,y1,'b-',t,y2,'g-',t,yex,'r--');            % plot t versus y1
    error1(hi)=y1(end)-yex(end);            % error at t=1.0 for Eulers
    error2(hi)=y2(end)-yex(end);            % error at t=1.0 for Imp. Euler
    idx=find(t==0.5);
    error3(hi)=y1(idx)-yex(idx);            % error at t=0.5 for Eulers
    error4(hi)=y2(idx)-yex(idx);            % error at t=0.5 for Imp. Euler

    disp(idx)
end
ratio1=error1(1:3)./error1(2:4);
ratio2=error2(1:3)./error2(2:4);

p1=log(ratio1)/log(2);
p2=log(ratio2)/log(2);

disp([p1' p2']);

ratio3=error3(1:3)./error3(2:4);
ratio4=error4(1:3)./error4(2:4);

p3=log(ratio3)/log(2);
p4=log(ratio4)/log(2);

disp([p3' p4']);