鉴于之间的依赖关系$z$和$t$，我们可以用时间导数重写第一个方程：

$$\partial_tU=\frac{c}{n}\left(\nabla_r^2U+\rho U\right)$$ . 现在我们定义$g=\left(\begin{array}{c}U\\\rho\end{array}\right)$. 这允许我们将方程组重写为： $$\partial_t g=\left(\begin{array}{cc}\frac{c}{n}\nabla_r^2&0\\0&0\end{array}\right)g+\frac{c}{n}\left(\begin{array}{c}\frac{1}{2}\\0\end{array}\right)g^T\left(\begin{array}{cc}0&1\\1&0\\\end{array}\right)g+\left(\begin{array}{c}0\\1\end{array}\right)a\left(g^T\left(\begin{array}{cc}1&0\\0&0\end{array}\right)g\right)^2$$

现在您可以为非线性方程组选择您最喜欢的时间步进算法：

$$\partial_tg=F(g,t,r)$$

例如，Crank Nicolson 将是一个可能的选择。现在没有必要忽略一些事情。请注意，您需要解决 Crank Nicolson 每一步的隐含部分的非线性。有关详细信息，请参阅此问题。

编辑

对于附加期限

$$\frac{c}{n} \beta|U|^2U$$ . 你需要添加 $$\frac{c}{n} \beta \left(\begin{array}{cc}1&0\\0&0\end{array}\right)g \left(g^T\left(\begin{array}{cc}1&0\\0&0\end{array}\right)g\right)$$ . $U$不应该在任何表达式中，我们要解决的领域$g$.

当然可以进一步简化表达式，但这应该足以作为一个起点。

如果您有更多条款$\partial_tU$, 简单地猜测$g$（这并不难），评估它，看看你是否达到了原始方程。

您正在做的是一种跳跃式方案，但任何时间步进方案都可以。但是由于您的问题是刚性和非线性的，您可能希望使用 B 稳定方法，以避免出于稳定性原因而不得不使用小时间步长。