What would be the best, or advised, method be to solve the following partial differential equation?
\begin{equation} i \hbar \frac{\partial}{\partial t} p_k = \Bigg(\epsilon_k^e + \epsilon_k^h - i \frac{\hbar}{T_2} \Bigg)p_k - (1- f_k^e-f_k^h) d_k E(t)+ieE(t) \cdot \nabla_k p_k \end{equation}
\begin{equation} \hbar \frac{\partial}{\partial t} f_k^{e} = -2 \textit{ Im}[d_kE(t)p_k^*] + eE(t) \cdot \nabla_k f_k^{e} \end{equation}
\begin{equation} \hbar \frac{\partial}{\partial t} f_k^{h} = -2 \textit{ Im}[d_kE(t)p_k^*] + eE(t) \cdot \nabla_k f_k^{h} \end{equation}
All quantities except $p_k(k,t)$, $f^e_k(k,t)$ and $f^h_k(k,t)$ are known. Furthermore, \begin{equation} \frac{\partial}{\partial t} f_k^e = \frac{\partial}{\partial t} f_k^h \end{equation} and \begin{equation} f_k^e = f_k^h \end{equation} have to be satisfied for all values of $k$. Initial conditions are all zero.
This closed set of coupled partial differential equations is known as the Semiconductor Bloch Equations. How would one go about numerically solving this? Is there a specific term for this type of equation within the field of (partial) differential equations?
MATLAB pdepe, seems to work for a slightly adapted version where the boundaries eventually approach zero. However, this adds up to a lot of computation time and is far from the ideal physical scenario which the equations should represent. It also seems to be the case that a quite large $k$ mesh resolution is required in order for the solver to be stable.
Any help is greatly appreciated.