%%  To solve a first order ODE:
syms y(t)
eqn = diff(y,t) == 1+2*t*y;
ySol(t) = dsolve(eqn);
disp(ySol(t));


%%  To solve a first order ODE initial value problem (IVP):
syms y(t) a
eqn = diff(y,t) == a*y;
cond = y(0) == 5;
ySol(t) = dsolve(eqn,cond);
disp(ySol(t));


%%  To solve a second order ODE:
syms y(t) a
eqn = diff(y,t,2) == a*y;
ySol(t) = dsolve(eqn);
disp(ySol(t));


%%  To solve a second order ODE initial value problem (IVP):
syms y(t) a b
eqn = diff(y,t,2) == a^2*y;
Dy = diff(y,t);
cond = [y(0)==b, Dy(0)==1];
ySol(t) = dsolve(eqn,cond);
disp(ySol(t));

%%  To solve a system of ODEs:
syms y(t) z(t)
eqns = [diff(y,t)==z, diff(z,t)==-y];
[ySol(t),zSol(t)] = dsolve(eqns);
disp([ySol(t),zSol(t)]);


%%  To solve a system of ODEs initial value problem (IVP):
syms y(t) z(t)
eqns = [diff(y,t)==z, diff(z,t)==-y];
cond = [y(0)==b, z(0)==1];
[ySol(t),zSol(t)] = dsolve(eqns,cond);
disp([ySol(t),zSol(t)]);



%%  To solve first order ODE numerically
function first_oder_ode
clc
% SOLVE  dt/dt = t - y^2  
% initial conditions: y(0) = 0

t=0:0.01:5;   % time scalex  <<<<  ************  CHANGE ME
initial_y=0;  <<<< ************  CHANGE ME

[t,y]=ode45( @rhs, t, initial_y);
disp([t,y])


plot(t,y);
xlabel('t'); ylabel('y');

    function dydt=rhs(t,y)
        dydt = t - y^2;   <<<< ************  CHANGE ME
    end
end


%%  To solve second order ODE numerically
clc
f = @(t,y) [y(2); t+y(2)-3*y(1)]; % define function f(t,y)
t0 = 0;  y0 = [1;-2];             % initial condition with vector y0
tf = 4;                            % final t-value
[ts,ys] = ode45(f,[t0,tf],y0);     % solve IVP

fprintf('y(t) = %g, y''(t) = %g\n',ys(end,1),ys(end,2))

disp('       t           y1          y2')
disp([ts,ys])                     % table of t and y values of solution

plot(ts,ys(:,1),'b')              % plot solution y(t)
title('Solution y(t) of IVP')
xlabel('t'); grid on


%%  eigenvalues:
A=[[3 4 ];[2 2]]
lambda = eig(A);
disp(lambda);

%%  eigenvectors:
A=[[3 4 ];[2 2]]
[V,D] = eig(A);
disp(V);

%%  Laplace Transforms:
syms t y
f = 1/sqrt(t);
disp(laplace(f));


syms a t
f = exp(-a*t);
L = laplace(f,y);
disp(L);

syms t s
syms a positive
disp(laplace(dirac(t-a),t,s));
disp(laplace(heaviside(t-a),t,s));


syms f(t) s
Df = diff(f(t),t);
laplace(Df,t,s)


%%  inverse Laplace Transforms:
syms t s
ilaplace((-1/3)*1/(s^2+4))


%  power series solution
syms y(x)

Dy = diff(y(x),x);

eqn = diff(y,2)+y==0;
dsolve(eqn,'ExpansionPoint',0,'Order',7)
 
eqn = diff(y,2)+y==x^2;
dsolve(eqn,'ExpansionPoint',0,'Order',10)




%%  HeavisideFunctions
syms x
plotAxes(-10,10,-10,20);
fplot((1/2)*x^2-heaviside(x-1)*((1/2)*(x-1)^2), [-1, 8],'b','LineWidth',5);