%% Numerical Simulation - Section 2.4, 2.5, 2.6
% One program to do it all: Euler, Improved Euler and Runge-Kutta
%
% Instructions: 
% (1) set the initial value (starting point of solution curve)
% (2) set the final time
% (3) set the step size h and compute the number of steps N,
% or alternatively set the number of steps N and compute h.
% (4) modify the function definition "f" to fit the DE
% (5) run the program (click the green triangle) to call the three
% numerical algorithms and to plot results. Label and save your graphs
% appropriately.

function [yE, yIE, yRK] = numerical()
% clear the command line
clc
% clear workspace memory
clear

% initial value
t0 = 0;
y0 = 0;

% final time
tf = 5;

% set step size h, then compute number of steps N
h = .1;
N = round((tf-t0)/h);

% or set number of steps N, then compute step size h
%N = 10;
%h = (tf-t0)/N;

% Call numerical method functions defined below
[tE,yE]   = Euler(@f,t0,y0,h,N);
[tIE,yIE] = Imp_Euler(@f,t0,y0,h,N);
[tRK,yRK] = RK(@f,t0,y0,h,N);

% Create a figure with plots of the outputs
figure(1)
plot(tE,yE,'r', tIE,yIE,'m', tRK,yRK,'b')
hold on
xlabel('t')
ylabel('y')
set(gca,'FontSize',12)
str = sprintf('Step Size h=%f',h)
title(str)
legend('Euler', 'Improved Euler', 'Runge-Kutta')

% Plot the exact solution, or just leave commented out.
%syms t

% Just use one of the plots below
%p1 = ezplot(sin(t), [t0 tf])
%p1 = ezplot(y0*exp(t), [t0 tf])
%p1 = ezplot(22*exp(t/5) - 5*t - 25,[t0 tf])

%set(p1, 'Color', 'k')
%hold off

end


%% Right hand side function of a first order ODE (dy/dt = f(t,y))
function [fOut] = f(t,y)
fOut = 10-0.2*y
%fOut = 10-0.2*y-0.27*y^(1.5)
%fOut = cos(t)
%fOut = y
%fOut = t+0.2*y
end

% function [err] = Error(ya,yb)
% err = abs((ya-yb)./yb)
% end


%% Euler
% Solve the IVP y'(t)=f(t,y) with y(t0)=y0. The function f is given as a
% symbolic function.
% Given: step size: h, and number of steps: N. 
function [t,y] = Euler(f,t0,y0,h,N)

t = zeros(1,N+1);
y = zeros(1,N+1);
t(1) = t0;
y(1) = y0;

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


%% Improved Euler
% Solve the IVP y'(t)=f(t,y) with y(y0)=t0. The function f is given as a
% symbolic function.
% Given: step size: h, and number of steps: N. 
function [t,y] = Imp_Euler(f,t0,y0,h,N)

t = zeros(1,N+1);
y = zeros(1,N+1);
t(1) = t0;
y(1) = y0;

for i=1:N                             
    k1 = f(t(i),y(i));
    k2 = f(t(i)+h,y(i)+h*k1);
    k = 0.5*(k1+k2);
    t(i+1) = t(i)+h;
    y(i+1) = y(i)+h*k;
end
end


%% Runge-Kutta
% Solve the IVP y'(t)=f(t,y) with y(t0)=y0. The function f is given as a
% symbolic function.
% Given: step size: h and number of steps: N. 
function [t,y] = RK(f,t0,y0,h,N)

t = zeros(1,N+1);
y = zeros(1,N+1);
t(1) = t0;
y(1) = y0;

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