Using MATLAB to solve a system of two PDEs encountered in Mathematical Biology I am in an undergraduate summer project and am studying a paper, I hope I could solve the following two PDEs with two boundary conditions for $y_1(t,a)$ and $y_2(t,a)$ using MATLAB.
$$\frac{\partial y_1}{\partial t}+\frac{\partial y_1}{\partial a}=f_1(t)y_1, \qquad t\geqslant0, \quad a\in[0,\tau_1],$$
with boundary condition
$y_1(t,0)=c_1$.
$$\frac{\partial y_2}{\partial t}+\frac{\partial y_2}{\partial a}=f_2(t)y_2, \qquad t\geqslant0,\quad a\in[0,\tau_2],$$
with boundary condition
$y_2(t,0)=c_2y_1(t,\tau_1)$.
Here, $y_1(t,a)$ and $y_2(t,a)$ are functions of $t$ and $a$. The functions $f_1(t)$ and $f_2(t)$ are known functions of $t$ but are rather complicated, and $\tau_1,\tau_2,c_1,c_2$ are known constants.
Since $y_1(t,a)$ is involved in the boundary condition of the second PDE for $y_2(t,a)$, I aim to solve the first PDE for $y_1(t,a)$ first, and then use the results in solving the second PDE for $y_2(t,a)$.
I am new in both PDEs and MATLAB, although I have spent several days searching for possible ways in doing it, I still have got no hints. It would be greatly appreciated if some hints or possible directions I could look into could be provided. Thanks in advance.
 A: First order PDEs can generally be solved by the method of characteristics. To solve your first PDE using this method, put
$$ \phi_z(x) = y_1(z + x, x)\ . $$
Then
\begin{eqnarray}
\frac{d\phi_z}{dx}(x)&=&\frac{\partial y_1}{\partial t}(z+x,x)+\frac{\partial y_1}{\partial a}(z+x,x)\\
&=& f_1(z+x)\,\phi_z(x)\ .
\end{eqnarray}
This is a first-order ODE for $\ \phi_z(x)\ $, which has the solution:
\begin{eqnarray}
\frac{\phi_z(x)}{\phi_z(0)}&=&e^{\int_0^xf_1(z +u)\,du}\\
&=& e^{\int_z^{x+z}f_1(u)\,du}\ .
\end{eqnarray}
Now, $\ \phi_z(0)=y_1(z, 0) = c_1\ $, and $\ y_1(t,a)=\phi_{t-a}(a)\ $, so
$$
y_1(t,a)=c_1 e^{\int_{t-a}^tf_1(u)\,du}\ .
$$
Your second PDE can now be solved similarly.  When solving the corresponding ODE, however, you will need to integrate over the interval from $\ x\ $ to $\ \tau_1\ $, rather than from $\ 0\ $ to $\ x\ $ as I did above.
A: The method of characteristics provides the solution in the domain $t\geq a \geq 0$. The characteristic curves along which the information $\frac{\text d}{\text d t} y_k = f_k y_k$ propagates are the straight parallel lines $x = t-t_0$ with $t_0 \geq 0$. Hence, we have
$$
y_1(t,a) = c_1 \exp\left(\int_{t-a}^t f_1(s)\, \text d s\right) = c_1 \exp\left(\int_{0}^a f_1(t-\tau)\, \text d \tau\right) .
$$
Similarly, we have for $t\geq a \geq 0$
$$
y_2(t,a) = c_2\, y_1(t,\tau_1) \exp\left(\int_{0}^a f_2(t-\tau)\, \text d \tau\right)
$$
where $y_1(t,\tau_1)$ with $t\geq \tau_1$ can be deduced from the previous step. Of course, this solution can then be implemented in Matlab, e.g. by using the Matlab function integral for the numerical integration of $f_k$. Alternatively, an iterative numerical resolution by finite-volume methods can be implemented.
