% Class 27
% More IVP ODE examples and adaptive methods and stiff problems
% 
close all
clc
% Example 1: Thanksgiving Dinner
yo = [18 10 10 2 5 5];
[tp,yp] = rk4sys(@SYS,[0 3],yo,0.01);
plot(tp,yp)
xlabel('Time (hr)')
ylabel('Quantity')
legend('T','S','M','G','P','J')
ylim([0 20])
% Use interpolation to solve for time value
tval = interp1(yp(:,4),tp,0)
% One more example of system of eqns.
% these can appear when you have a second order problem
[tp,yp] = ode45(@SYS2,[0 2],[pi/2 0])
plot(tp,yp(:,1))
% A few examples of STIFF problems
% 23.8 from the book
dydt = @(t,y) 5*(y-t^2);
fan = @(t) t^2 + 0.4*t +0.08;
fplot(fan,[0 5])
[tp,yp] = rk4sys(dydt,[0 5],0.08,0.03125)
ylim([0 30])
hold on
plot(tp,yp)
% Now try ode45
[tp,yp] = ode45(dydt,[0 5],0.08);
plot(tp,yp)
% Now try ode15s
[tp,yp] = ode15s(dydt,[0 5],0.08);
plot(tp,yp)
[tp,yp] = rk4sys(dydt,[0 5],0.08,0.0001);
plot(tp,yp,'o')
legend('Analytical','RK 0.031 step','ode45','ode15s','RK 0.0001 step')
% Fogler example, of system of ODE in reaction engineering
tspan = [0 0.5];
y0 = [0.021 0.0105 0];
[tp, yp] = ode45(@SYS3,tspan,y0);
figure
plot(tp,yp)
legend('Hydrogen','Mesitylene','Xylene')
xlabel('Tau (hr)')
ylabel('C (lbmol/ft^3)')
function dy = SYS3(t,y)
% t = tau (residence time)
% y(1) = hydrogen
% y(2) = m_input
% y(3) = xylene
k1 = 55.2;
k2 = 30.2;
dy = [-k1*y(1)^(.5)*y(2)-k2*y(3)*y(1)^(.5);...
    -k1*y(2)*y(1)^(.5);....
    k1*y(2)*y(1)^(.5)-k2*y(3)*y(1)^0.5];
end

function dy = SYS2(t,y)
% t = time
% y(1) = theta
% y(2) = h
dy = [y(2);...
    -(9.8)/.6*sin(y(1))];
end


function dy = SYS(t,y)
% t = time (independent variable)
% y(1) = T
% y(2) = S
% y(3) = M
% y(4) = G
% y(5) = P
% y(6) = J
dy = [-1-y(2)/2*t;...
    -1-y(1)/5*t;...
    -3-y(4)/2*t;...
    -1-y(3)/5*t;...
    -20/y(1)*t;...
    (y(1)+y(2)+y(3)+y(4)+y(5))/10*t];
end