%% Modifying Euler's Method
% First we recall Euler's Method for the example \frac{dy}{dt}=t+y, y(0)=1. We seek an approximation to y(1).
%

clear
f = @(t,y) t+y;         % RHS of equation
a = 0;   b = 1;         % Interval [a,b]
N = 10;                  % Number of subintervals
h = (b-a)/N;            % Stepsize
t = linspace(a,b,N+1);
y(1) = 1;               % Initial value

for k = 2:N+1
    y(k) = y(k-1) + h*f(t(k-1),y(k-1));
end

%
% The exact solution is y(t)=2e^t-t-1 and y(1)=2e-2=3.4366... .
%

plot(t,y,'or')

%

sol = @(x) 2*exp(x)- x - 1;
hold on
fplot(sol,[0,1],'b')
hold off
Euler_err=abs(sol(1)-y(N+1))

%
%% Modified Euler

N = 10;                 % Number of subintervals
h = (b-a)/N;            % Stepsize
t = linspace(a,b,N+1);
y(1) = 1;               % Initial value

for k=1:N
    k1    = y(k) +  h*f(t(k),y(k));
    y(k+1)= y(k) + h/2*(f(t(k),y(k)) + f(t(k)+h,k1));
end
ModifiedEuler_err=abs(sol(1)-y(N+1))

%
%% Huen Method

N = 10;                 % Number of subintervals
h = (b-a)/N;            % Stepsize
t = linspace(a,b,N+1);
y(1) = 1;               % Initial value

for k=1:N
    k1    = y(k) +  2/3*h*f(t(k),y(k));
    y(k+1)= y(k) + h/4*(f(t(k),y(k)) + 3*f(t(k)+2/3*h,k1));
end
Huen_err=abs(sol(1)-y(N+1))

%
%% General Second Order Methods
% Using Taylor series expansions and y''(x)=f_x(x,y)+f(x,y)f_y(x,y), y(x+h)=y(x)+hy'(x)+\frac{h^2}{2}y''(x)+\cdots \\
% y(x+h)=y(x)+h \left[ b_1f(x,y) + b_2f(x+c_2h,y+a_{2,1}h f(x,y)) \right]
% Comparisons give b_1+b_2=1,\quad b_2c_2=\frac 12, \quad b_2a_{2,1}=\frac 12.
% Let a_{2,1}=\alpha.Then, c_2=a_{2,1}=\frac 1{2b_2}, \quad b_2=\frac 1{2\alpha}, \quad b_1=1-\frac 1{2\alpha},
% Butler Tableau: \left( \begin{array}{ccc} 0 \vert & 0 & 0 \\c_2 \vert & a_{2,1} & 0  \\ \hline\quad \vert & b_1 & b_2 \end{array} \right) = \left( \begin{array}{ccc} 0 \vert & 0 & 0 \\ \alpha \vert & \alpha & 0  \\ \hline\quad \vert & 1-\frac{1}2{\alpha} & \frac{1}2{\alpha}  \end{array} \right).
%
%% Example of Runge-Kutta Schemes:
% y_{n+1}=y_n+h(b_1k_1+b_2k_2+\cdots+b_mk_m),
% k_1=f(x_n,y_n),
% k_2=f(x_n+c_2h,y_n+a_{21}hk_1), ...
% k_m=f(x_n+c_mh,y_n+a_{m,1}hk_1+\cdots a_{m,m-1}hk_{m-1}),
% Examples:
% \left( \begin{array}{cccc} Method &\alpha & b_1 & b_2  \\ \hline
% Huen & 1 & \frac 12 & \frac 12 \\
% Midpoint& \frac 12 & 0 & 1 \\
% Ralston \; I & \frac 34 & \frac 13 & \frac 23 \\
% Ralston \; II & \frac 23 & \frac 14 & \frac 34
% \end{array} \right)
%  Note: There are some differences as to which is Huen and which is Ralston. Read about the controversy making our Ralston II the correct one as the error satisfies |E|<\left(4\left|\frac 16 - \frac{\alpha}4\right| + \frac 13\right)ML^2. This is minimized for \alpha=\frac 23.

alpha = 2/3;   beta = alpha;  b2 = 1/2/alpha; b1 = 1-b2; % 
N = 10;                 % Number of subintervals
h = (b-a)/N;            % Stepsize
t = linspace(a,b,N+1);
y=zeros(1,N+1);
y(1) = 1;               % Initial value

for k=1:N
    k1 = f(t(k),y(k));
    k2 = f(t(k)+alpha*h,y(k)+beta*k1*h);
    y(k+1)= y(k) + h*(b1*k1+b2*k2);
end
SecondOrder_err=abs(sol(1)-y(N+1))

%
% It would be nice to plot |E|/ML^2 vs \alpha and point out the special cases above.
%

M=10;
alpha=zeros(1,M);
E=zeros(1,M);
for j=1:M
    alpha(j)=j/M;
    beta = alpha(j);  b2 = 1/(2*alpha(j)); b1 = 1-b2; % 
    N = 10;                 % Number of subintervals
    h = (b-a)/N;            % Stepsize
    t = linspace(a,b,N+1);
    y(1) = 1;               % Initial value
    
    for k=1:N
        k1 = f(t(k),y(k));
        k2 = f(t(k)+alpha(j)*h,y(k)+beta*k1*h);
        y(k+1)= y(k) + h*(b1*k1+b2*k2);
    end
    E(j) = abs(sol(1)-y(N+1));
end
plot(alpha,E,'o')

%
%% Fourth Order Runge-Kutta

N = 50;                 % Number of subintervals
h = (b-a)/N;            % Stepsize
t = linspace(a,b,N+1);
y(1) = 1;               % Initial value

for i=1:N
    k1=h*f(t(i),y(i));
    k2=h*f(t(i)+h/2,y(i)+k1/2);
    k3=h*f(t(i)+h/2,y(i)+k2/2);
    k4=h*f(t(i)+h,y(i)+k3);
    y(i+1)=y(i)+(k1+2*k2+2*k3+k4)/6;
end
RKerr=abs(sol(1)-y(N+1))