clear all

% Matlab project 2: solving initial value problems using Euler, modified
%                   Euler and Runge-Kutta 4th order methods
%   solves IVP dy/dx=x+y, y(0)=1

%N=10;
x0=0; y0=1; xf=5;
%h=(xf-x0)/N;
h=0.01;
N=(xf-x0)/h;

[Xe,Ye]=eulerODE(x0,y0,xf,N); % use Euler method

[Ximpe,Yimpe]=impeuler(x0,y0,xf,N); % improved Euler method

[Xrk,Yrk]=rk(x0,y0,xf,N);   % Runge-Kutta 4th order method

xex(1)=x0; yex(1)=y0;
for i=1:N
   xex(i+1)=x0+i*h;
   yex(i+1)=exact_sol(xex(i+1)); 
end


% Compute errors between exact and computed numerically solutions

for i=1:N+1
   erre(i)=abs(yex(i)-Ye(i)); 
   errimpe(i)=abs(yex(i)-Yimpe(i)); 
   errrk(i)=abs(yex(i)-Yrk(i)); 
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

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

hold on

plot(Xe,Ye,'b','Linewidth',2,'Marker','o')

plot(Ximpe,Yimpe,'r-.','Linewidth',2,'Marker','*')

plot(Xrk,Yrk,'g','Linewidth',2,'Marker','s')

plot(xex,yex,'k:','Linewidth',2,'Marker','.')
 
legend('Euler','Improved Euler','Runge-Kutta','exact solution',2)

xlabel('x')
ylabel('y')

title('Solution to dy/dx=x+y, y(0)=1')

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

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

hold on

plot(xex,erre,'b','Linewidth',2,'Marker','o')
plot(xex,errimpe,'r','Linewidth',2,'Marker','*')
plot(xex,errrk,'g','Linewidth',2,'Marker','s')
title('Absolute error')
legend('Euler','Improved Euler','RK')
xlabel('x')
ylabel('Absolute error')

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

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

hold on

plot(xex,erre./abs(yex),'b','Linewidth',2,'Marker','o')
plot(xex,errimpe./abs(yex),'r','Linewidth',2,'Marker','*')
plot(xex,errrk./abs(yex),'g','Linewidth',2,'Marker','s')

xlabel('x')
ylabel('Relative error')

title('Relative error')
legend('Euler','Improved Euler','RK',2)