%% Solve a BVP by means of Symbolic Math Toolbox
% Let us consider the following PDE problem, i.e. the 1D wave equation with initial and boundary conditions:
% \left\{\begin{array}{lll}
% \frac{\partial^2u}{\partial t^2} = \frac{\partial^2u}{\partial x^2}, &  0<x<L, & t > 0 \\
% u(x,0^+)   = \frac{x\sin(3x)}{6} \\
% \frac{\partial u}{\partial t}(x,0^+) = \frac{\sin(3x)}{6} + \frac{x\cos(3x)}{2} \\
% u(0,t) = \frac{t\sin(3t)}{6} \\
% u(L,t) = \frac{(L+t)\sin[3(L+t)]}{6}
% \end{array}\right.
% Its analytical solution is the function: u(t,x)=\frac{(x+t)\sin[3(x+t)]}{6}.
% In order to solve the problem, we apply the Laplace Transform method to the PDE problem with respect to the variable t. It converts the PDE problem to an ODE problem, and produces the following BVP (Boundary Value Problem), where the unknown function is the Laplace Transform of the solution of the PDE: U(x,s)={\cal L}_t[u(t,x)].
% \left\{\begin{array}{lcll}
% U'' &=& s^2\,U - (s\,x+1)\frac{\sin(3x)}{6} - x\frac{\cos(3x)}{2},  & 0 < x < L \\[8pt]
% U(0,s) &=& \frac{s}{(s^2 + 9)^2}, \\[8pt]
% U(L,s) &=& \frac{s\cos(3L)-3\sin(3L)}{(s^2 + 9)^2} + \frac{(sL+1)\sin(3L)+3L\cos(3L)}{6(s^2 + 9)}
% \end{array}\right.
% The complex variable s, introduced by the Laplace Transform, is considered to be a parameter and will be ignored in solving the BVP problem symbolically.
% The analytical solution of the BVP problem is: U(x,s) = \frac{(s\,x+1)\sin(3x) + 3x\cos(3x)}{6(s^2+9)} + \frac{s\cos(3x) - 3\sin(3x)}{(s^2+9)^2}.

%% wave_BVP.mlx

clear; clc; close all

%
%% Analytical Laplace Transform of u(t,x): U(x,s) = \frac{(s\,x+1)\sin(3x) + 3x\cos(3x)}{6(s^2+9)} + \frac{s\cos(3x) - 3\sin(3x)}{(s^2+9)^2}

%{
solution U(x,s) for x in]0,L[, s Complex
of the ODE BVP: diff(U,x,2) == s^2*U - (s*x+1)*sin(3*x)/6 - x*cos(3*x)/2
conditions:     U(0,s) == s/(s^2+9)^2
                U(L,s) == (s*cos(3*L)-3*sin(3*L))/(s^2+9)^2 + ((sL+1)*sin(3*L)+3*L*cos(3*L))/(6*(s^2+9))
%}
syms x L real
syms s
syms U(x) % s is considered as a parameter

ODE = diff(U,x,2) == s^2*U - (s*x+1)*sin(3*x)/6 - x*cos(3*x)/2;

cond1 = U(0) == s/(s^2+9)^2;
cond2 = U(L) == (s*cos(3*L)-3*sin(3*L))/(s^2+9)^2 + ((s*L+1)*sin(3*L)+3*L*cos(3*L))/(6*(s^2+9));
conds = [cond1, cond2]; % boundary conditions

%% solve the BVP symbolically
Usol = dsolve(ODE,conds)

% known analytical solution
Utrue = ((s*x+1)*sin(3*x)+3*x*cos(3*x))/(6*(s^2+9)) + (s*cos(3*x)-3*sin(3*x))/(s^2+9)^2;
% comparison between solutions
fprintf('\nCheck the symbolic solution is equal to the analytical one: Usol - Utrue = '); disp(simplify(Usol - Utrue))

%

% compute the Laplace Transform U(x,s) for s=1
Usol1 = simplify(subs(Usol,s,1))
% plot the Laplace Transform U(x,s) for s=1
figure(1); clf
fplot(Usol1,[-3*pi,3*pi])
xlabel('$x$','Interpreter','latex','FontSize',14)
title('Laplace Transform $U(x,s)={\cal L}_t\left[u(x,t)\right],$ for $s=1$','Interpreter','latex','FontSize',16)

%
%% Analytical solution of the PDE problem: u(t,x)=\frac{(x+t)\sin[3(x+t)]}{6}
% Once the Laplace Transform U(x,s) has been computed by solving the BVP, in order to find symbolically the solution of the PDE problem, we need to apply the Inverse Laplace Transform u(t,x)={\cal L}^{\,-1\,}[U(x,s)] by means of the symbolic function ilaplace():

syms t real
u = simplify(ilaplace(Usol,s,t)) 
% known analytical solution
uTRUE = (x+t)*sin(3*(x+t))/6;
% comparison between solutions
fprintf("\nCheck the symbolic Inverse Laplace Transform u(t,x) is equal to the analytical one: u - uTRUE = "); disp(simplify(u - uTRUE))

%

%% display the surface u(t,x)
Tmax=10; L=1;
figure(2); clf
fsurf(u,[0 8 0 1]) % 0 <= t <= 8,  0 <= x <= 1
xlabel('$t$','Interpreter','latex','FontSize',14)
ylabel('$x$','Interpreter','latex','FontSize',14)
title('Inverse Laplace Transform: $u(t,x)={\cal L}^{\,-1}_t\left[U(x,s)\right]$','Interpreter','latex','FontSize',16)
%% create the animated wave
Nt=50; Nx=20;
[tj,xi]=meshgrid(linspace(0,Tmax,Nt),linspace(0,L,Nx));
uij=double(subs(u,{t,x},{tj,xi}));
umax=max(max(double(uij)));
umin=min(min(double(uij)));

figure(5),clf
for j=1:Nt
    plot(xi(:,1),uij(:,j))
    axis([0 L umin umax])
    xlabel('$x$','Interpreter','latex','FontSize',14)
    title(['Wave $u(t,x)$ for $t=$' num2str(tj(1,j))],'Interpreter','latex','FontSize',14)
    pause(0.25)
end

%
%% Insert an animated gif
%