clear all;
close all;

% Euler Method

t0=0;
u0=1;

h = input(' Enter the value of the time step ')
n=input(' Enter the number of time steps ')

t=t0;
u=u0;

T(1)=t0;
sol(1)=u0;

for i=1:n
    k1=f(t,u);
    u=u+h*k1;
    t=t+h;
T(i+1)=t;
sol(i+1)=u;
end;

% Improved Euler Method

t=t0;
u=u0;

for i=1:n
    k1=f(t,u);
    k2=f(t+h,u+h*k1);
    u=u+0.5*h*(k1+k2);
    t=t+h;
    sol_iEm(i+1)=u;
end;

sol_iEm(1)=u0;


% Runge-Kutta Method

t=t0;
u=u0;

for i=1:n
    k1=f(t,u);
    k2=f(t+0.5*h,u+0.5*h*k1);
    k3=f(t+0.5*h,u+0.5*h*k2);
    k4=f(t+h, u+h*k3);
    u=u+(h/6)*(k1+2*k2+2*k3+k4);
    t=t+h;
    sol_RK(i+1)=u;
end;

sol_RK(1)=u0;

plot(T,sol_exact(T),'r',T,sol,T,sol_iEm,T,sol_RK,'k');
title('Numerical methods for ODEs')
legend('exact', 'Euler', 'Modified Euler', 'RK','Location','NorthEast')


% The error of these methods

for i=1:n+1
    E_Euler(i)=abs(sol(i)-sol_exact(T(i)));
    E_iEm(i)=abs(sol_iEm(i)-sol_exact(T(i)));
    E_RK(i)=abs(sol_RK(i)-sol_exact(T(i)));
end;

figure
plot(T,E_Euler, T,E_iEm, T, E_RK);
title('Error functions');
legend('Euler', 'IEM', 'RK');


%axis ([0 T(n+1) 0 sol(n+1)],'square');