%% Laplace transform for ODE with piecewise forcing function

clearvars
syms s t Y  

% We want to solve the initial value problem
%         y''+3y'+2y=f(t),\qquad y(0)=0,\quad y'(0)=0
% where the forcing function is defined as a piecewise function:
%          f(t) = \cases{ t & for $0\le t<1$ \cr
% 0& for $t\ge 1$
% }
%% Step 1: Rewrite f(t) using step function u(t)
% We can rewrite a piecewise defined function using the unit step function (Heaviside function) u(t) as follows
% f(t) = \cases{ f_1(t) & for $0\le t<c$ \cr
% f_2(t) & for $t\ge c$
% }  \\
% f(t) = f_1(t)+ u(t-c)\bigl(f_2(t)-f_1(t)\bigr)
% Here we obtain
%         f(t) = t+u(t-1)(0-t)

f = t + heaviside(t-1)*(-t);
fplot(f,[0 6]); title('forcing function f(t)')

%% Step 2: Find Laplace transform F=\mathcal{L}[f]
% We have
%         F(s) = \mathcal{L}\left[ t
% \right]  +
% \mathcal{L}\left[u(t-1)\,(-t)
% \right]
% For the first term we get \mathcal{L}\left[t\right]=\frac1{s^2}
% For the second term we use the following rule:
%%          \mathcal{L}\left[ u(t-c)g(t)
% \right] = e^{-cs}\mathcal{L}\left[g(t+c)
% \right]
% Here g(t)=-t, so g(t+1) = -(t+1).
% Now we find \mathcal{L}\left[ g(t+1)
% \right] =
% \mathcal{L}\left[ -1-t
% \right] = -\frac1s-\frac1{s^2}.
% This gives
%         F(s)=\frac1{s^2}+e^{-s}\left(
% -\frac1s-\frac1{s^2}
% \right)

F = laplace(f)

%% Step 3: Find Laplace transform of left hand side and solve for Y(s)
% We have
%          Y_1(s)=\mathcal{L}[y']=s\cdot Y(s)-y(0)=s\cdot Y(s) \\
% Y_2(s)=\mathcal{L}[y'']=s\cdot Y_(s)-y'(0)=s^2\cdot Y(s)

Y1 = s*Y-0;
Y2 = s*Y1-0;

% We obtain
%         (s^2+3s+2)Y(s)=F(s) \\
% Y(s) = \frac1{(s^2+3s+2)}
% \left[\frac1{s^2}+e^{-s}\left(
% -\frac1s-\frac1{s^2}
% \right)
% \right]

Sol = simplify( solve(Y2+3*Y1+2*Y==F,Y) ,1000)
collect(Sol,'exp')

%% Step 4: Find inverse Laplace transform
% We obtained Y(t) = R_1(s)+e^{-s}R_2(s) with rational functions
%          R_1(s) =  \frac1{(s^2+3s+2)}\cdot
% \frac1{s^2}    \\
% R_2(s) = \frac1{(s^2+3s+2)}
% \left(
% -\frac1s-\frac1{s^2}
% \right)
% = \frac{-1}{(s+2)(s+1)}\cdot\frac{s+1}{s^2}
% =\frac{-1}{s+2}\cdot\frac1{s^2}
% We need to find the partial fraction decomposition for R_1 and R_2:

R1 = partfrac(1/(s^2*(s + 1)*(s + 2)))
R2 = partfrac((-1/(s^2*(s + 2))))

% The solution y(t) is given by
%         y(t) = \mathcal{L}^{-1}\left[Y(t)\right]
% = \mathcal{L}^{-1} \left[R_1(s)
% \right]
% +\mathcal{L}^{-1} \left[e^{-cs}R_2(s)
% \right]
% It is clear how to find \mathcal{L}^{-1} \left[R_1(s)
% \right]:

ilaplace(R1)

% For \mathcal{L}^{-1} \left[e^{-cs}R_2(s)
% \right] we use the rule
%%         \mathcal{L}\left[
% u(t-c)g(t-c)
% \right] = e^{-cs}G(s)
% We first find g(t)=\mathcal{L}^{-1}\left[ R_2(s)
% \right]

g = ilaplace(R2)

% We then obtain
%        \mathcal{L}^{-1}\left[
% e^{-s}G(s)
% \right] = u(t-1)g(t-1) = u(t-1)\left( \frac14-\frac14e^{-2(t-1)}-\frac{(t-1)}2
% \right)
% In Matlab we can just apply ilaplace to Sol:

sol = ilaplace(Sol)
fplot(sol,[0 6]); title('solution y(t)'); axis auto

% We can collect the terms with 'heaviside':

sol = collect(sol,'heaviside')

% We can rewrite this as a "piecewise defined function":

sol = rewrite(sol,'piecewise')

% We can see that the solution for t>1 converges to zero as t\to\infty.