I would like to obtain analytic solutions to the following PDE system: \begin{equation} \rho_t + D(\lambda)\,\rho_\lambda = A(\lambda) \rho, \tag{1} \end{equation} with $\rho = (\rho_0,\rho_1)^T$, $D$ is the diagonal matrix \begin{equation} D(\lambda) = \begin{pmatrix} -(1+a)\lambda & 0 \\ 0 & p_b-\lambda\end{pmatrix}, \end{equation} and the matrix $A$ depends on $\lambda$ alone, through component functions of the form $ \frac{\lambda^2 + a \lambda + b}{c \lambda + d}$.

I am aware of the fact that questions of this type have been asked before on MathOverflow. However, system $(1)$ has certain features which make it, in my opinion, worthwhile to devote a new question to it.

- System $(1)$ is non-homogeneous (i.e. $A \neq 0$), so it cannot be interpreted as a conservation equation.
- The left hand side of $(1)$ is in diagonal form, while $A$ has no nonzero components. In other words, $D$ and $A$ are not simultaneously diagonalizable.

While I am aware of the fact that 'There is no generally applicable method of characterstics for first order systems' (@Igor Khavkine's comment on this question), I would hope that the characteristic surface spanned by the two characteristics defined by the left hand side of $(1)$ could be used to solve $(1)$ by a generalisation of the method of characteristics.

Since system $(1)$ does not have constant coefficients, previous questions such as this, this, and this do not apply, as far as I can see. Moreover, the nonvanishing right hand side of $(1)$ makes the situation qualitatively different from this and this question. Also, system $(1)$ cannot be reduced in a manner demonstrated in this question or this question.

Any ideas would be highly appreciated!