function ODEpractice
% Coded by NFR on 11.16.2017
% This code shows some examples of ODEs
%
% Problem 22.1
%
D = @(t,y) y.*t.^3-1.5.*y;
% Part a) Solve analytically
Yfunc = @(t) exp(t.^4/4-1.5.*t);
%fplot(Yfunc,[0 2])
tvec = linspace(0,2,100);
yvec = Yfunc(tvec);
plot(tvec,yvec)
hold on
xlabel('t')
ylabel('y')
m = 10;
% Euler's Method (using Chapra's code)
[t,y] = eulode(D,[0 2],1,0.5)
[t2,y2] = eulode(D,[0 2],1,0.25)
plot(t,y,t2,y2)

% Solve with midpoint formula (2 pt Rk, huen method)
ti = 0;
tf = 2;
h = 0.5;
t = (ti:h:tf)'; 
n = length(t);
y0 = 1;
y = y0*ones(n,1); %preallocate y to improve efficiency
for i = 1:n-1 %implement Euler's method
    yeuler = y(i) + D(t(i),y(i))*(t(i+1)-t(i)); % Euler method
    yhalf = (yeuler + y(i))/2;
    thalf = (t(i+1)+t(i))/2;
    slopehalf = D(thalf,yhalf);
    y(i+1) = y(i) + slopehalf*(t(i+1)-t(i));
end
plot(t,y)
% Let's use 4 point Runga Kutta Method with h = 0.5
h = 0.5;
[tp,yp] = rk4sys(D,[0 2],1,0.5);
plot(tp,yp)
hold off
legend('Analytical','Euler w/ 0.5 step size','Euler w/ 0.25 step size',...
    'Mid-point formula with h = 0.5', '4pt RK with h = 0.5')
% Now, let's look at systems of ODEs
% Problem 22.7 from the text
ti = 0;
tf = 0.4;
h = 0.1;
y0 = 2;
z0 = 4;
DY = @(t,y) -2.*y+5.*exp(-t);
DZ = @(t,y,z) -y.*z.^2/2;
tvec = ti:h:tf;
ns = length(tvec);
% Preallocate for outputs
Yvec = y0*ones(ns,1);
Zvec = z0*ones(ns,1);
for i = 1:ns-1
    Yvec(i+1) = Yvec(i)+DY(tvec(i),Yvec(i))*h;
    Zvec(i+1) = Zvec(i)+DZ(tvec(i),Yvec(i),Zvec(i))*h;
end
plot(tvec,Yvec,tvec,Zvec)

% Looks good, but let's check with the RK algorithm
tspan = [0 .4];
y0 = [2;4];
h = 0.1;
[tp,yp] = rk4sys(@SYS,tspan,y0,h);
hold on
plot(tp,yp(:,1),tp,yp(:,2))
hold off
legend('Y','Z','Yrk','Zrk')



end

function dy = SYS(t,y)
% y(1) = Y
% y(2) = Z
dy = [-2.*y(1)+5.*exp(-t);...
    -y(1).*y(2).^2/2];
end