%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Solves IVP y'=x+y subject to IC y(x0)=y0, x0=0, y0=1
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clear all

x0=0; xf=5;
y0=1;
%N=200;
%h=(xf-x0)/N;
h=0.1;
N=(xf-x0)/h;


X=zeros(N+1,1);
Y=zeros(N+1,1);
yex=zeros(N+1,1);
error=zeros(N+1,1);
error_rel=zeros(N+1,1);

xold=x0; yold=y0;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   This solution does not save intermediate results
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for i=1:N
   xnew=xold+h;
   ynew=yold+h*(xold+yold);
   disp(['i=',num2str(i),', x=',num2str(xnew),', y=',num2str(ynew)])
   xold=xnew; yold=ynew;
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   Compute exact solution at x points
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

X(1)=x0; Y(1)=y0;   % initial condition
yex(1)=-X(1)-1+2*exp(X(1)); % same as y0
error(1)=0;
error_rel(1)=0;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   This part saves all intermediate points for plotting
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

for i=1:N
  X(i+1)=X(i)+h;
  Y(i+1)=Y(i)+h*(X(i)+Y(i));

  yex(i+1)=-X(i+1)-1+2*exp(X(i+1));
  error(i+1)=abs(Y(i+1)-yex(i+1));  % absolute error
  error_rel(i+1)=error(i+1)/yex(i+1);   % relatiive error
end


figure(1); clf(1)
set(gca,'FontSize',12); set(gca,'box','on')
plot(X,yex,'r','Linewidth',2)
hold on
plot(X,yex,'rx','Linewidth',2)
plot(X,Y,'b-.','Linewidth',2)
plot(X,Y,'bo','Linewidth',2)
title('Numerical solution')
legend('exact','','Euler','')

figure(2); clf(2)
set(gca,'FontSize',12); set(gca,'box','on')

plot(X,error,X,error,'k','Linewidth',2)
title('Absolute error')

figure(3); clf(3)
set(gca,'FontSize',12); set(gca,'box','on')

plot(X,error_rel,X,error_rel,'k','Linewidth',2)
title('Relative error')