%% SOLVING DIFFERENTIAL EQUATIONS


%% [Exercise 1] system solution for t=3:
clear all
close all 
clc

t=3;
initial_conditions=[0,1,1];
solution=model_1(t,initial_conditions);

%% [Exercise 2]:system solution for several time intervals
%Notes:An example of a nonstiff system is the system of equations describing
%the motion of a rigid body without external forces.
clear all
close all 
clc

%time_interval=[0 12];
time_interval=0:12;
initial_condition=[0 1 1];
%options = odeset('RelTol',1e-4,'AbsTol',[1e-4 1e-4 1e-5]);
%[T,Y] = ode45(@model_1,time_interval,initial_condition,options);

[T,Y] = ode45(@model_1,time_interval,initial_condition);

plot(T,Y(:,1),'-',T,Y(:,2),'-.',T,Y(:,3),'.')

%% [Exercise 3]: simple example
clear all
close all 
clc

%Equation: y'=5y-3

time_interval=0:0.1:5;
initial_condition=[1];


[T,Y] = ode45(@model_4,time_interval,initial_condition);

subplot(1,2,1);plot(T,Y(:,1),'--*b');
xlabel('Time')
ylabel('Y [numerical]')

%Exact Solution:
c=initial_condition;
t=time_interval;
Y_exact=(1/5)*exp(5*t+5*c)+(3/5);



subplot(1,2,2);plot(t,Y_exact,'--*r')
xlabel('Time')
ylabel('Y [exact]')
%% [Exercise 4]:system solution for several time intervals
%Notes:An example of a stiff system is provided by the van der Pol equations
%in relaxation oscillation. The limit cycle has portions where the solution
%components change slowly and the problem is quite stiff, alternating with regions
%of very sharp change where it is not stiff.
clear all
close all 
clc

time_interval=[0 3000];
initial_condition=[2 0];

%[T,Y] = ode45(@model_2,time_interval,initial_condition);
[T,Y] = ode15s(@model_2,time_interval,initial_condition);

subplot(1,2,1);plot(T,Y(:,1),'-ob')
subplot(1,2,2);plot(T,Y(:,2),'-*r')




%% Solving Second order Differential equations in matlab:
%Consider an 80Kg paratrooper falling from 600 meters.
%The trooper is accelerated by gravity, but decelerated by drug on the
%parachute.
%Equation: m(d2y/dt2)=-mg+4/5 V^2
clear all
close all
clc

time_interval=linspace(0,15,100);%[0 15]; %seconds
initial_conditions=[600 0];

[T,Y] = ode45(@model_5,time_interval,initial_conditions);

figure(1)
subplot(1,2,1);plot(T,-Y(:,1),'-*b')
ylabel('height(m)')
xlabel('time(s)')
subplot(1,2,2);plot(T,-Y(:,2),'-*r')
ylabel('velocity(m)')
xlabel('time(s)')


%% FUNCTIONS USED:
% Note: To use the functions you need to cut them from here and paste them
% in a new function file separately !!!
% %------------------------------------------------------------
% function [ dy ] = model_1( t,y )
% 
% dy=zeros(3,1);
% dy(1)=y(2)*y(3);
% dy(2)=-y(1)*y(3);
% dy(3)=-0.51*y(1)*y(2);
% 
% end
% %------------------------------------------------------------
% function [ dy] = model_2( t,y )
% 
% dy = zeros(2,1);    % a column vector
% dy(1) = y(2);
% dy(2) = 1000*(1 - y(1)^2)*y(2) - y(1);
% 
% end
% %------------------------------------------------------------
% function [ dy ] = model_4( t,y )
% 
% dy=zeros(1,1);
% dy=5*y(1)-3;
% 
% end
% %------------------------------------------------------------
% function rk=model_5(t,ic)
% %ic:initial conditions (2x1)
% %ic(1):height
% %ic(2);velocity
% 
% 
% mass=80;
% g=9.81;
% 
% rk=[ic(1);-g+(4/5)*ic(2).^2/mass];
% 
% end
% %------------------------------------------------------------
% 
% function dx=model_5(t,y)
% 
% mass=80;
% g=9.81;
% 
% dy(1)=y(1)
% dy(2)=-g+(4/5)*y(2).^2/mass;
% 
% end