clc;clear;close all
% Model.
% dVdt = qin-c*V-s*Tg-(1-s)*Tl;
% dQdT = 1/Q;
% Tg^3 + Q*Tg^2+Tg-1.2 = 0; s*(Tl^3 + Q*Tl^2+Tl-1.2);(this can be written without s*, but the result is the same) for example 11.
% T^3 + Q*T^2+T-1.2 = 0; for example 12.


%% Example 11. Solves the ODE system using ode15s while the algebraic are solved internally with root function
clear 
clc
y0 = [1 1.5];
tspan = [0 10];
options = odeset('Stats','on');
[t,y,stats] = ode15s(@odefun11,tspan,y0,options);
figure
plot(t,y,'Linewidth',1)
legend('x','y','Location','northwest')
%axis([0 5 0 3.5])
xlabel('Time','Fontsize',20,'Fontname','Timesnewroman')
ylabel('Differential variables','Fontsize',20,'Fontname','Timesnewroman')
set(gca,'fontsize',20)
grid on

%% Example 12. Solves the DAE system using ode15s 
% at t=0, uses an imaginary root for Tl, since ode15s computes only one
% root which is closest to the initial condition
% This leads to imaginary values of V however. 

y0 = [1 1.5 0.5583 -1.0292+1.044i]; 
%y0 = [1 1.5 0.5583 0.5583]; 
%y0 = [1 1.5 0.5583 1i]; 


tspan = [0 10];
M = [1 0 0 0; 0 1 0 0; 0 0 0 0; 0 0 0 0];
options=odeset('Mass',M);
[t,y] = ode15s(@odefun12,tspan,y0,options);
figure
plot(t,y(:,1:2),'Linewidth',1)
legend('x','y','Location','northwest')
%axis([0 5 0 3.5])
xlabel('Time','Fontsize',20,'Fontname','Timesnewroman')
ylabel('Differential variables','Fontsize',20,'Fontname','Timesnewroman')
set(gca,'fontsize',20)
grid on
