% Solve du/dt=u-2t/u fpr t>0 with u(0)=1
% Exact solution is u=sqrt(2t+1)
% We use explicit Euler, improved Euler, modified Euler and RK4
% We use a range of time-step sizes to show convergence rates

% right-hand-side function
f=@(t,u) u-2*t./u;

% exact solution
u=@(t) sqrt(2*t+1);

% solution at t=1
u1=u(1);

Rmax=11;
% set up vectors to store errors
e1=zeros(Rmax,1); e2=zeros(Rmax,1); e3=zeros(Rmax,1); e4=zeros(Rmax,1); 
% set up vector of N values
Ns=2.^((1:Rmax)+2);

% run schemes with each value of N in turn
for r=1:Rmax
    N=Ns(r);
    U=explicit_euler(f,N); e1(r)=abs(U(end)-u1);
    U=improved_euler(f,N); e2(r)=abs(U(end)-u1);
    U=modified_euler(f,N); e3(r)=abs(U(end)-u1);
    U=myRK4(f,N);          e4(r)=abs(U(end)-u1);
end

loglog(Ns,e1,'bx-',Ns,e2,'rx-',Ns,e3,'gx-',Ns,e4,'kx-',Ns,1./Ns,'b--',Ns,1./Ns.^2,'r--',Ns,1./Ns.^4,'k--')
xlabel('N')
ylabel('errors')
legend('Explicit Euler','Improved Euler','Modified Euler','RK4','O(1/N)','O(1/N^2)','O(1/N^4)','location','SouthWest')

f=@(t,u) -25*u;
U=myRK4(f,9);


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% function for explicit Euler
function U=explicit_euler(f,N)
t=linspace(0,1,N+1);
dt=t(2)-t(1);
U=zeros(1,N+1); U(1)=1;
for n=1:N
    U(n+1)=U(n)+dt*f(t(n),U(n));
end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% function for improved Euler
function U=improved_euler(f,N)
t=linspace(0,1,N+1);
dt=t(2)-t(1);
U=zeros(1,N+1); U(1)=1;
for n=1:N
    k1=f(t(n),U(n));
    k2=f(t(n)+dt,U(n)+dt*k1);
    U(n+1)=U(n)+dt*(k1+k2)/2;
end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% function for modified Euler
function U=modified_euler(f,N)
t=linspace(0,1,N+1);
dt=t(2)-t(1);
U=zeros(1,N+1); U(1)=1;
for n=1:N
    k1=f(t(n),U(n));
    k2=f(t(n)+dt/2,U(n)+dt*k1/2);
    U(n+1)=U(n)+dt*k2;
end
end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% function for RK4
function U=myRK4(f,N)
t=linspace(0,1,N+1);
dt=t(2)-t(1);
U=zeros(1,N+1); U(1)=1;
for n=1:N
    k1=f(t(n),U(n));
    k2=f(t(n)+dt/2,U(n)+dt*k1/2);
    k3=f(t(n)+dt/2,U(n)+dt*k2/2);
    k4=f(t(n)+dt,U(n)+dt*k3);
    U(n+1)=U(n)+dt*(k1+2*k2+2*k3+k4)/6;
end
end