%% Solving initial value problems with the Laplace transform
% Let
%           Y(s)=\mathcal{L}[y], \qquad
% Y_1(s)=\mathcal{L}[y'], \qquad
% Y_2(s)=\mathcal{L}[y'']
% Since \mathcal{L}[f'] = sF(s)-f(0) we have
%           Y_1(s)=sY(s)-y(0),\qquad
% Y_2(s) = sY_1(s)-y'(0)

clearvars
syms s t Y       % define symbolic variables s,t,Y

%% Example 1: y'-3y=e^{2t},\quad y(0)=4
% Take the Laplace transform of the right-hand side function f(t)=e^{2t}:

f = exp(2*t);
F = laplace(f)

% Find Y_1 = sY(s)-y(0):

Y1 = s*Y-4;

% Solve the equation Y_1(s)-3Y(s)=F(s) for Y(s):

Sol = solve(Y1-3*Y==F,Y)

% Find the partial fraction decomposition by hand: Y(s) = \frac{A}{s-2}+\frac{B}{s-3}, multiply by (s-2)(s-3):
%        1+4(s-2) = A(s-3)+B(s-2)
% Plugging in s=2 gives 1=A(-1) or A=-1.
% Plugging in s=3 gives 5 = B\cdot 1 or B=5.
% Find the partial fraction decomposition using Matlab:

partfrac(Sol)

% Now we can apply the inverse Laplace transform:  Using \mathcal{L}^{-1}[\frac1{s-2}]=e^{2t},\quad
% \mathcal{L}^{-1}[\frac1{s-3}]=e^{3t} we obtain
%         y(t) = 5e^{3t}-e^{2t}
% Use ilaplace to find the inverse Laplace transform:

sol = ilaplace(Sol)

%
%% Example 2: y''+3y'+2y=e^{-t},\quad y(0)=4,\quad y'(0)=5
% Here Y(s)=\frac{p(s)}{q(s)} and q(s) has a double root -1.

f = exp(-t);
F = laplace(f)

% Find Y_1 = sY(s)-y(0):

Y1 = s*Y-4;

% Find Y_2 = sY_1(s)-y'(0):

Y2 = s*Y1-5;

% Solve the equation Y_2(s)+3Y_1(s)+2Y(s)=F(s) for Y(s):

Sol = solve(Y2+3*Y1+2*Y==F,Y)

% Find the partial fraction decomposition by hand:
% We first need the roots of the characteristic polynomial: p(s)=s^2+3s+2=0 gives s=-1 and s=-2, so we have s^2+3s+2=(s+1)(s+2). The denominator function for Y(s) will be
%      (s+1)(s^2+3s+2)=(s+1)^2(s+2)
% Hence we want to find a partial fraction decomposition
%       Y(s) = \frac{A}{s+1}+\frac{B}{(s+1)^2}+\frac{C}{s+2}
% Multiplying by (s+1)^2(s+2) gives
%       (4s+17)(s+1)+1 = A(s+1)(s+2) + B(s+2) + C(s+1)^2        (*)
% Plugging in s=-1 gives 1 = B(-1+2), so B=1.
% Plugging in s=-2 gives 9\cdot (-1)+1=C(-2+1)^2, so C=-8.
% It remains to find A.
% Method 1: Take the derivative of (*):
%        4(s+1)+(4s+17)1 = A(s+1+s+2)+B+C\cdot 2(s+1)
% In the derivative plug in s=-1. This gives  13= A\cdot 1+B, so A=13-B=13-1=12.
% Method 2: In (*), plug in a different value, e.g., s=0: this gives 18=A\cdot 2+B\cdot 2+ C\cdot 1, A = (18-2B-C)/2=(18-2+8)/2=12
% Find the partial fraction decomposition using Matlab:

partfrac(Sol)

% Now we can apply the inverse Laplace transform:  Using \mathcal{L}^{-1}[\frac1{s+1}]=e^{-t},\quad
% \mathcal{L}^{-1}[\frac1{(s+1)^2}]=te^{-t},\quad
% \mathcal{L}^{-1}[\frac1{s+2}]=e^{-2t} we obtain
%          y(t) = 12 e^{-t}+te^{-t}-8e^{-2t}
% Use ilaplace to find the inverse Laplace transform:

sol = ilaplace(Sol)

%% Example 3: y''-2y'+5y=3,\quad y(0)=1,\quad y'(0)=1
% Here Y(s)=\frac{p(s)}{q(s)} and q(s) has complex roots 1+2i,\; 1-2i.
% So far we just used F=laplace(f)
% The full form of the command is F=laplace(f,t,s) where t is the variable for f, and s is the variable for F.

f = 3;
F = laplace(f,t,s)    % since there is no explicit t in f we need to use laplace(...,t,s)
Y1 = s*Y-1;
Y2 = s*Y1-0;

% Here is Y(s):

Sol = solve(Y2-2*Y1+5*Y==F,Y)

% Taking the inverse Laplace transform gives the solution y(t):

sol = ilaplace(Sol)

% If we want to find this by hand we first need the partial fraction decomposition. There are two ways to do this:
%% Method 1: Partial fraction decomposition (real version)

partfrac(Sol)

% We have \mathcal{L}[\frac1s]=1, so \mathcal{L}[\frac{3/5}s]=\frac35.
% We have s^2-2s+5=(s-1)^2+2^2=(s-a)^2+b^2 since the roots are a\pm ib. We want to write the term as
%        \frac{A(s-a)+B}{(s-a)^2+b^2}
% yielding
%        \mathcal{L}
% \left[
% \frac{A\cdot(s-a)+\frac{B}b\cdot b}{(s-a)^2+b^2}
% \right] = A e^{at}\cos(bt)+ \frac{B}be^{at}\sin(bt)
% Here \frac25s -\frac45= \frac25(s-1)+\frac25 - \frac45= \frac25(s-1)-\frac25, so
%          \mathcal{L}
% \left[
% \frac{\frac25(s-1)-\frac15\cdot 2}{(s-1)^2+2^2}
% \right]
% =
% \frac25e^t\cos(2t)-\frac15\sin(2t)
%% Method 2: Partial fraction decomposition (complex version)

partfrac(Sol,'FactorMode','full')

% We get complex roots a+ib,\; a-ib. If these occur once (as in this example) the partial fraction decomposition has terms
%         \frac{C}{s-(a+ib)}+\frac{\overline{C}}{s-(a-ib)}
% with complex C=A+iB, \overline{C}=A-iB.
% After we multiply by the denominators we can find the complex number C by plugging in s=a+ib.
% Then we obtain
%          v := \mathcal{L}^{-1}\left[
% \frac{C}{s-(a+ib)}
% \right]
% =C e^{(a+ib)t} = (A+iB)e^{at}\bigl(\cos(bt)+i\cdot \sin(bt)\bigr)
%         \mathcal{L}^{-1}\left[
% \frac{C}{s-(a+ib)}+\frac{\overline{C}}{s-(a-ib)}
% \right] = v + \overline{v} = 2\text{Re}\,v = 2e^{at}\bigl(B\cos(bt)-A\sin(bt)\bigr)
% If the complex roots a+ib,\; a-ib occur twice (as in Example 4 below) we will get additional terms
%         \frac{D}{(s-(a+ib))^2}+\frac{\overline{D}}{(s-(a-ib))^2}
% with complex D=E+iF,\; \overline{D}=E-iF.
% We obtain
%         v := \mathcal{L}^{-1}\left[
% \frac{D}{(s-(a+ib))^2}
% \right]
% =D\, t\,e^{(a+ib)t} = (E+iF)t\,e^{at}\bigl(\cos(bt)+i\cdot \sin(bt)\bigr)
%         \mathcal{L}^{-1}\left[
% \frac{D}{(s-(a+ib))^2}+\frac{\overline{D}}{(s-(a-ib))^2}
% \right] = v + \overline{v} = 2\text{Re}\,v = 2\,t\,e^{at}\bigl(E\cos(bt)-F\sin(bt)\bigr)
% In example 3 we have simple roots a+ib,\; a-ib and C = \frac15+ i\frac1{10}, so
%          \mathcal{L}\left[
% \frac{\frac15+i\frac1{10}}{s-1-2i}+\frac{\frac15-i\frac1{10}}{s-1+2i}
% \right] = \frac25e^t\cos(2t)
% -\frac15e^t\sin(2t)
%
%% Example 4: y''+y=\sin(t)+\cos(t),\quad y(0)=1,\quad y'(0)=2
% Here Y(s)=\frac{p(s)}{q(s)} and q(s) has double complex roots i,\;-i.

f = sin(t)+cos(t);
F = laplace(f)
Y1 = s*Y-1;
Y2 = s*Y1-2

% We obtain
%        (s^2Y-s-2)+Y= \frac{s+1}{s^2+1}
%        Y = \frac{s+2 + \frac{s+1}{s^2+1}}{s^2+1}

Sol = simplify( solve(Y2+Y==F,Y) )

% Partial fraction decomposition by hand: (complex version)
%        Y = \frac{s+2 + \frac{s+1}{s^2+1}}{s^2+1}
% = \frac{C}{s-i}+\frac{\overline{C}}{s+i}+\frac{D}{(s-i)^2} +
% \frac{\overline{D}}{(s+i)^2}
% multiply by (s^2+1)^2 = (s-i)^2(s+i)^2:
%        (s^2+1)(s+2)+(s+1) = C(s+i)(s^2+1) + \overline{C}(s-i)(s^2+1)
% + (s+i)^2 D + (s-i)^2 \overline{D}
% plug in s=i:
%        0+(i+1) = (i+i)^2 D gives i+1=-4D, so D=-\frac14-\frac14i
% take the derivative:
%       (s^2+1)+2s(s+2)+1 = C\bigl((s^2+1)+(s+i)2s\bigr)
% + \overline{C}\bigl((s^2+1)+(s-i)2s\bigr)
% +2(s+i)D + 2(s-i)\overline{D}
% plug in s=i:
%        2i(i+2)+1 = C2i\cdot2i+2\cdot 2i\cdot D
%        -2+4i+1 = -4C-4i(-\frac14-\frac14i)
%         -1+4i=-4C-i+1,\qquad
% 4C = -5i+2, \qquad C = \frac12-\frac54i
% Partial fraction decomposition with Matlab: (complex version)

partfrac(Sol,'FactorMode','full')
sol = ilaplace(Sol)