## 5 Numerical Solution to the Inverse Problem  
When the four boundary conditions are chosen, we have  
\[u(x,y,t)=U(x,y,t)+Z(x,y,t), \tag{20}\]  
where \(U(x,y,t)\) is the solution to problem (18), (19) when \(F\equiv 0\) and \(Z(x,y,t)\)--when \(F\
eq 0\) but all boundary conditions and initial condition are also zero. Then, \(Z(x,y,t)\) takes the form  
\[Z(x,y,t)=\\ \int_{0}^{X}\int_{0}^{Y}\int_{0}^{t}Z(x,y;\alpha,\gamma,t-\tau)\exp \biggl{(}-\beta\tau-\mu x-\eta\sqrt{\frac{a}{c}}y\biggr{)}f(\tau)\delta\biggl{(} \alpha-x_{0},\sqrt{\frac{a}{c}}(\gamma-y_{0})\biggr{)}\mathrm{d}\alpha\mathrm{ d}\gamma\mathrm{d}\tau\\ =\int_{0}^{t}Z(x,y,x_{0},y_{0},t-\tau)\exp\biggl{(}-\beta\tau-\mu x -\eta\sqrt{\frac{a}{c}}y\biggr{)}f(\tau)\mathrm{d}\tau. \tag{21}\]  
Now, when (4), (5) (respectively, (1), (2)) is a problem with unknown time-dependent strength \(f(t)\), it is referred to as an inverse problem. To reconstruct the unknown function in (21), the overposed data condition (7) is applied:  
\[u(x_{1},y_{1},t)=U(x_{1},y_{1},t)+Z_{1}(x_{1},y_{1},t)=Y(t)\exp\biggl{(}-\beta t -\mu x_{1}-\eta\sqrt{\frac{a}{c}}y_{1}\biggr{)},\ \ 0\leq t\leq T.\]  
Thus, we obtain a Volterra integral equation for \(f(t)\):  
\[g(t)=Y(t)\exp\biggl{(}-\beta t-\mu x_{1}-\eta\sqrt{\frac{a}{c}} y_{1}\biggr{)}-U(x_{1},y_{1},t)\\ =\int_{0}^{t}Z(x_{1},y_{1};x_{0},y_{0},t-\tau)\exp\biggl{(}-\beta \tau-\mu x_{1}-\eta\sqrt{\frac{a}{c}}y_{1}\biggr{)}f(\tau)\mathrm{d}\tau,\ \ 0\leq t\leq T. \tag{22}\]  
Since the function \(g(t)\) is known, the problem (22) has a unique solution \(f(t)\in\mathbb{L}^{2}(0,T)\); see, e.g., [22]. We use a collocation numerical method in order to derive the solution to the first-kind Volterra Equation (22), with the convolution kernel  
\[R(x_{1},y_{1};x_{0},y_{0},t-\tau)=Z(x_{1},y_{1};x_{0},y_{0},t-\tau).\]  
We introduce a time uniform mesh \(t_{k}=k\tau,\ k=0,1,\ldots,M,\ M\tau=T\). Then, we seek the approximate solution  
\[f(t)\approx f_{m}(t)=\sum_{m=1}^{M}f_{m}\phi_{m}(t), \tag{23}\]  
where \(\phi_{m}(t)\) is the \(m\)-th base function defined by  
\[\phi_{m}(t)=\left\{ \begin{array}{ll}\dfrac{t-t_{m-1}}{t_{m}-t_{m-1}},&t_{m-1 }\leq t\leq t_{m},\\ \dfrac{t_{m+1}-t}{t_{m+1}-t_{m}},&t_{m}\leq t\leq t_{m+1},\\ 0,&\text{elsewhere}.\end{array} \right.\]  
We note that \(\{\phi_{m}(t)\}_{m=1}^{M}\) is an orthonormal set in \(\mathcal{C}[0,T]\). By placing (23) in (22), at successive time \(t=t_{k},k=1,\ldots,M\), we obtain  
\[g(t_{k})=\sum_{m=1}^{k}f_{m}\int_{0}^{t_{k}}R(x_{1},y_{1};x_{0},y_{0},t_{k}- \tau)\phi_{m}(\tau)\mathrm{d}\tau=\sum_{m=1}^{k-1}f_{m}a_{k-m+1}+f_{m}a_{1},\ k=1, \ldots,M,\]  
where  
\[a_{1}=\int_{0}^{t_{1}}R(x_{1},y_{1};x_{0},y_{0},t_{1}-\tau)\phi_{1}(\tau) \mathrm{d}\tau,\]  
\[a_{k}=\int_{0}^{t_{2}}R(x_{1},y_{1};x_{0},y_{0},t_{k}-\tau)\phi_{1}(\tau) \mathrm{d}\tau,\quad k=2,\ldots,M,\]  
and the respective integrals are calculated with a suitable numerical procedure [23]. Now, let us consider the linear algebraic system of equations  
\[Af=g,\ \ f=(f_{1},\ldots,f_{k})^{\top},\ \ g=(g_{1},\ldots,g_{k})^{\top}, \tag{24}\]  
which is obtained by (23) and (24), such that \(A=\mathds{R}^{M\times M}\) is a lower-triangular Toeplitz matrix given by  
\[A=\left( \begin{array}{cccc}a_{1}&0&\cdots&0\\ a_{2}&a_{1}&\cdots&0\\ \vdots&\vdots&\ddots&\vdots\\ a_{M}&a_{M-1}&\cdots&a_{1}\end{array} \right),\ \ a_{k}>0,\ \ k=1,\ldots,M.\]  
In order to solve (24), we implement the Tikhonov regularization method, and for this, we define the discrete functional  
\[J(f)=\sum_{m=1}^{M}\Bigg{\{}\big{(}\sum_{k=1}^{m}(a_{m-k+1}f_{k}-g_{k})\big{)} ^{2}+\alpha\big{(}\sum_{k=1}^{m}l_{m-k+1}f_{k}\big{)}^{2}\Bigg{\}}, \tag{25}\]  
where \(\alpha>0\) is a predefined regularization parameter and \(L\) is a matrix of type \(A\) instead of identity one \(I\) in the sequential Tikhonov regularization algorithm. To estimate the unknown parameter \(f\) using the least-squares method, we minimize the function \(J(f)\) as shown in (25).