% This program implements Euler's method, Heun's method, and Runge-Kutta
% method. It analyses the approximation error as function of step size h.

% Problem: See Example 1 in sections 8.1 , 8.2 , 8.3 , 
% equations (5), (6), (7).
%      y'=1-t+4y , y(0)=1
% Exact solution: Phi(t) = t/4 - 3/16 + 19/16 exp(4t)

% Compare three methods for the following step sizes:
% Euler's method: hEuler
% Heun's method: hHeun = 2*hEuler
% Runge-Kutta: hRK = 4*hEuler
% These are chosen to have the same number of f-computations over the
% interval [0,1]

clear all;
tmax = 1;

hEulervector = [0.000001,0.000002,0.000004,0.000005,0.00001,0.00002,...
    0.00004,0.00005,0.0001,0.0002,0.0004,0.0005,0.001,0.002,0.005,...
    0.01,0.05];


imax = length(hEulervector);
%Allocate memory
hh = zeros(imax,1);
yEuler = zeros(imax,1);
yHeun = zeros(imax,1);
yRK = zeros(imax,1);
yexact = zeros(imax,1);


for indk = 1:length(hEulervector);
    
hEuler = hEulervector(indk);
hHeun = 2*hEuler;
hRK = 4*hEuler;

hh(indk) = hEuler;

% Euler's method:
h = hEuler;
y = 1;
for t=0:h:tmax-h
    f = 1-t+4*y;
    y = y+h*f;
end
yEuler(indk) = y;

% Heun's method:
h = hHeun;
y = 1;
for t=0:h:tmax-h
    f1 = 1-t+4*y;
    f2 = 1-(t+h)+4*(y+h*f1);
    y = y + 0.5*(f1+f2)*h;
end
yHeun(indk) = y;

% Runge-Kutta
h = hRK;
y=1;
for t=0:h:tmax-h
    k1 = 1-t+4*y;
    k2 = 1-(t+h/2)+4*(y+h*k1/2);
    k3 = 1-(t+h/2)+4*(y+h*k2/2);
    k4 = 1-(t+h)+4*(y+h*k3);
    y = y + (k1+2*k2+2*k3+k4)*h/6;
end
yRK(indk) = y;

% Exact solution
% Phi(t) = t/4 - 3/16 + 19/16 exp(4t)
t = tmax;
yexact(indk) = t/4-3/16 + 19*exp(4*t)/16;

end


% Analysis of Results
figure(1)
plot(hh,yEuler-yexact,'r',hh,yHeun-yexact,'g',hh,yRK-yexact,'b')
xlabel('h'),ylabel('Error'),legend('Euler','Heun','Runge-Kutta','Location','NorthWest')
%print -djpeg NumericalAnalysisComp.jpg

figure(2)
loglog(hh,abs(yEuler-yexact),'r',hh,abs(yHeun-yexact),'g',hh,abs(yRK-yexact),'b')
xlabel('h'),ylabel('Error'),legend('Euler','Heun','Runge-Kutta','Location','NorthWest')
%print -djpeg NumericalAnalysisCompLogLogScale.jpg
