Similar to the question I asked on the Math Stackexchange (in this question), I am interested in solving a nonlinearly coupled system of PDEs.
The system of PDEs look like this: $$\begin{align*}\partial_t p_{11} &= E_p p_{13} - E_p p_{31} \\ \partial_t p_{12} &= E_p p_{13} - E_p p_{32} - p_{12} (\Delta_{31}-\Delta_{32}) \\ \partial_t p_{13} &= E_p p_{11} + E_c p_{12} - E_p p_{33} + p_{13} \Delta_{31} \\ \partial_t p_{21} &= E_p p_{23} - E_c p_{31} + p_{21}(-\Delta_{31} + \Delta_{32}) \\ \partial_t p_{22} &= E_c p_{23} - E_c p_{32} \\ \partial_t p_{23} &= E_p p_{21} + E_c p_{22} - E_c p_{33} + p_{23} \Delta_{31} + p_{23} (-\Delta_{31} + \Delta_{32}) \\ \partial_t p_{31} &= - E_p p_{11} - E_c p_{21} + E_p p_{33} - p_{31} \Delta_{31} \\ \partial_t p_{32} &= - E_p p_{21} - E_c p_{22} - E_c p_{33} + p_{23} \Delta_{31} + p_{23} (-\Delta_{31} + \Delta_{32}) \\ \partial_t p_{33} &= E_p p_{13} - E_c p_{23} + E_p p_{31} + E_c p_{32} \end{align*}$$
$$ \frac{\partial E_p}{\partial z } + \frac{1}{c} \frac{\partial E_p}{\partial t } = i k p_{13}(t, z) $$
Each of these p's ($p_{11}, p_{12}, ..., p_{33}$) are a function of t and z. Additionally $E_p$ is a function of t and z. And $\Delta_{31, 32, 33}$ are just numerical constants.
I'm trying with the boundary conditions: $$\text{p11}(0,z)=1,\\ \text{p12}(0,z)=\text{p13}(0,z)=\text{p21}(0,z)= \text{p22}(0,z)=\text{p23}(0,z)=\text{p31}(0,z)= \text{p32}(0,z)= \text{p33}(0,z)=0,\\ E_p(t,0)=e^{-(t-t_0)^2}- e^{-(t_0)^2}\\ E_p(0,z)=0 $$
Is there a way that I can solve this numerically in Mathematica?
When I try the following code:
NDSolve[{
Derivative[1, 0][p11][t, z] == -(1/2) I Ep[t, z] (p13[t, z] - p31[t, z]),
Derivative[1, 0][p12][t, z] ==
1/2 I (-10 p13[t, z] + Ep[t, z] p32[t, z]),
Derivative[1, 0][p13][t, z] == -(1/2) p13[t, z] -
1/2 I (10 p12[t, z] + Ep[t, z] (p11[t, z] - p33[t, z])),
Derivative[1, 0][p21][t, z] == -(1/2) p21[t, z] -
1/2 I (Ep[t, z] p23[t, z] - 10 p31[t, z]),
Derivative[1, 0][p22][t, z] == -5 I (p23[t, z] - p32[t, z]),
Derivative[1, 0][p23][t, z] == -(1/2)
I (Ep[t, z] p21[t, z] + 10 p22[t, z] - I p23[t, z] -
10 p33[t, z]),
Derivative[1, 0][p31][t, z] ==
1/2 I (10 p21[t, z] + I p31[t, z] +
Ep[t, z] (p11[t, z] - p33[t, z])),
Derivative[1, 0][p32][t, z] ==
1/2 I (Ep[t, z] p12[t, z] + 10 p22[t, z] + I p32[t, z] -
10 p33[t, z]),
Derivative[1, 0][p33][t, z] ==
1/2 I (10 p23[t, z] + Ep[t, z] (p13[t, z] - p31[t, z]) -
10 p32[t, z] + 2 I p33[t, z]),
Derivative[0, 1][Ep][t, z] +
Derivative[1, 0][Ep][t, z] == 1/2 I p12[t, z], p11[0, z] == 1,
p12[0, z] == 0, p13[0, z] == 0, p21[0, z] == 0, p22[0, z] == 0,
p23[0, z] == 0, p31[0, z] == 0, p32[0, z] == 0, p33[0, z] == 0,
Ep[t, 0] == E^-t^2}, {ρ11[t, z], ρ12[t, z], ρ13[t,
z], ρ21[t, z], ρ22[t, z], ρ23[t, z], ρ31[t,
z], ρ32[t, z], ρ33[t, z], Ep[t, z]}, {z, 0, 1}, {t, 0,
1}]
Which looks like this:
Mathematica returns the error message: NDSolve::femnonlinear: Nonlinear coefficients are not supported in this version of NDSolve.
Any ideas?
EDIT: I modified the boundary conditions for $E_p$ as suggested in the comments. (Previously I only had E_p(t,0)=e^{-(t-t_0)^2}, and now I have added an initial condition and modified the previous B.c. to make it consistent.).