% Improved Euler method
%
%  solves eqn  2.5.6  y'=-2*x*y, y(0)=2 
%  the exact answer  y(x) = 2*exp(-x^2)
%  
clear all;
f=inline('-2*x*y'); yx=inline('2*exp(-x.^2)');
x0=0.0; y0=2.0; 
n=10;  xx=zeros(1,n+1); yy=zeros(1,n+1); 
xn=1.0; 
h=(xn-x0)/n;
xx(1)=x0; yy(1)=y0;
for i=1:n
    x=xx(i); y=yy(i); k1=f(x,y);
    x=x+h; u=y+h*k1;  % predictor step
    k2=f(x,u); k=0.5*(k1+k2);
    y=y+h*k;          % corrector step 
    xx(i+1)=x; yy(i+1)=y;
end

%figure; plot(xx,yy,'b*'), 
%title('Improved Euler method: Approximate solution');
x=linspace(x0,xn,100); 
figure; plot(x,yx(x),'r-',xx,yy,'b*'); 
title('Improved Euler method: Exact and computed solutions')