% This program implements Euler's method, Heun's method, and Runge-Kutta
% method.

% 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: h = 0.05
% Heun's method: h = 0.1
% Runge-Kutta: h = 0.2
% These are chosen to have the same number of f-computations over the
% interval [0,1]

tmax = 1;
hEuler = 0.05;
hHeun = 2*hEuler;
hRK = 4*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 = 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 = 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 = y;

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

fprintf(1,'At time t = %f: \n',tmax);
fprintf(1,'yEuler (h = %f) = %f \n',hEuler,yEuler);
fprintf(1,'yHeun (h = %f) = %f \n',hHeun,yHeun);
fprintf(1,'yRunge-Kutta (h = %f) = %f \n',hRK,yRK);
fprintf(1,'Exact Solution = %f\n',yexact);
