% Class 25
% Coded by Nigel F. Reuel om 11.19.19
%
% Book problem 22.1 
%
% Analytical solution
F_an = @(t) exp(t^4/4-1.5*t);
fplot(F_an,[0 2])
hold on
% Numerical methods to solve this ODE
% Euler's method BIG STEP
h =0.5;
yo = 1;
tsteps = 0:h:2;
nsteps = length(tsteps);
dydt = @(t,y) y*t^3-1.5*y;
% Preallocate a vector to store the solution
Yval = zeros(nsteps,1);
% Store our initial point in this solution vector
Yval(1,1) = yo;
for i = 1:nsteps-1
    t = tsteps(i);
    slope = dydt(t,yo);
    ynew = yo + h*slope;
    yo = ynew;
    Yval(i+1,1) = ynew;
end
% Compare our euler solution to the analytical
plot(tsteps,Yval,'r.','MarkerSize',20)
%
% Euler's method SMALL STEP
h =0.25;
yo = 1;
tsteps = 0:h:2;
nsteps = length(tsteps);
dydt = @(t,y) y*t^3-1.5*y;
% Preallocate a vector to store the solution
Yval = zeros(nsteps,1);
% Store our initial point in this solution vector
Yval(1,1) = yo;
for i = 1:nsteps-1
    t = tsteps(i);
    slope = dydt(t,yo);
    ynew = yo + h*slope;
    yo = ynew;
    Yval(i+1,1) = ynew;
end
% Compare our euler solution to the analytical
plot(tsteps,Yval,'k.','MarkerSize',20)
%
% As an example of all the other RK methods
% we will do the midpoint form
%
% Part c
%
h =0.5;
yo = 1;
tsteps = 0:h:2;
nsteps = length(tsteps);
dydt = @(t,y) y*t^3-1.5*y;
% Preallocate a vector to store the solution
Yval = zeros(nsteps,1);
% Store our initial point in this solution vector
Yval(1,1) = yo;
for i = 1:nsteps-1
    t = tsteps(i);
    % This is a two point method to find
    % the increment function (phi) 
    % aka better estimate of the slope
    k1 = dydt(t,yo);
    k2 = dydt(t+h,yo+k1*h);
    slope = 1/2*k1 + 1/2*k2; % this is phi
    ynew = yo + h*slope;
    yo = ynew;
    Yval(i+1,1) = ynew;
end
% Compare our euler solution to the analytical
plot(tsteps,Yval,'g.','MarkerSize',20)
%
% Part (d) 4 pt method RK
%
h = 0.5;
yo = 1;
tsteps = 0:h:2;
nsteps = length(tsteps);
dydt = @(t,y) y*t^3-1.5*y;
% Preallocate a vector to store the solution
Yval = zeros(nsteps,1);
% Store our initial point in this solution vector
Yval(1,1) = yo;
for i = 1:nsteps-1
    t = tsteps(i);
    % This is a four point method to find
    % the increment function (phi) 
    % aka better estimate of the slope
    k1 = dydt(t,yo);
    k2 = dydt(t+1/2*h,yo+1/2*k1*h);
    k3 = dydt(t+1/2*h,yo+1/2*k2*h);
    k4 = dydt(t+h,yo+k3*h);
    slope = 1/6*(k1+2*k2+2*k3+k4); % this is phi
    ynew = yo + h*slope;
    yo = ynew;
    Yval(i+1,1) = ynew;
end
% Compare our euler solution to the analytical
plot(tsteps,Yval,'m.','MarkerSize',20)
% Obviously the 4 pt method is superior
% there must be a built in equation, indeed!
%
% CHAPRA has one - rk4sys
% MATLAB has ode45
%
[tp,yp] = rk4sys(dydt,[0 2],1,0.5)
plot(tp,yp,'ko','MarkerSize',10)
%
[tp,yp] = ode45(dydt,[0 2],1)
plot(tp,yp,'bo','MarkerSize',10)
%
% Solve systems of ODEs
% Using a built in function, it is easiest
% to code both functions as a local (sub) function
[tp,yp] = ode45(@SYS,[0 0.4],[2 4])
hold off
plot(tp,yp(:,1),tp,yp(:,2))
legend('y values','z values')

function dydt = SYS(t,y)
% y(1) = y
% y(2) = z
dydt = [-2*y(1)+5*exp(-t);
    -y(1)*y(2)^2/2];
end