%  Euler's method to solve IVP
%      dy/dx = f(x,y)
%      y(x_0) = y_0
%  from x=xbegin to x=xend.
%  f(x,y) and the exact solution for comparison are to be entered in files
%  frhs.m and exact_solu.m

xbegin=input('enter the value of xbegin: ');
xend=input('enter the value of xend: ');
n=input('enter the number of steps: ');

h=(xend-xbegin)/n;

x = linspace(xbegin,xend,n+1);
y = zeros(n+1,1);   % approximate solution from Euler's

%  specify the initial condition:
y(1) = exact_solu(x(1));

yimp = y;           % approximate solution from improved Euler
yexact = y;         % exact solution

for k=1:n
    y(k+1) = y(k) + h*frhs(x(k),y(k));
    k1 = frhs(x(k),yimp(k));
    u = yimp(k) + h*k1;
    k2 = frhs(x(k+1),u);
    yimp(k+1) = yimp(k) + 0.5*h*(k1+k2);
    yexact(k+1) = exact_solu(x(k+1));
end

fprintf('\n*** Euler method printout ***\n\n');

for k=1:n+1
fprintf('x = %6.4f  y_exact = %6.4f  y_app = %6.4f\n',x(k),yexact(k),y(k));
end

% print scaled errors:
fprintf('\nEuler error (scaled by h) at x=xend: %6.4f\n',abs(yexact(n)-y(n))/h);
fprintf('Improved Euler error (scaled by h^2) at x=xend: %6.4f\n',abs(yexact(n)-yimp(n))/h^2);

hold on

plot(x,yexact,'-k');
xlabel('x'); ylabel('y'); title('IVP solutions (exact and approximate)');
%axis([xbegin xend -6 10]);

pause

plot(x,y,'xr');

pause

plot(x,yimp,'+g');

legend('exact','Euler','improved');
