I would like to numerically solve the following system of coupled nonlinear differential equations:
$$ -\frac{\hbar^2}{2m_a} \frac{\partial^2}{\partial x^2}\psi_a + V_{ext}\psi_a + \left( g_a |\psi_a|^2 + g_{ab} |\psi_b|^2 \right)\psi_a = \mu_a \psi_a $$
$$ -\frac{\hbar^2}{2m_b} \frac{\partial^2}{\partial x^2}\psi_b + V_{ext}\psi_b + \left( g_b |\psi_b|^2 + g_{ab} |\psi_a|^2 \right)\psi_b=\mu_b\psi_b $$
where $\hbar$, $m_a$, $m_b$, $g_a$, $g_b$, $g_{ab}$ are known coefficients and $V_{ext}$ is a known function of $x$, i.e.: $$ V_{ext}= -P \left[\cos\left(\frac{3}{2}\, \frac{x}{L}\, 2\pi \right)\right]^2 $$ The unknowns are eigenfunctions $\psi_a(x)$, $\psi_b(x)$ and the eigenvalues $\mu_a$ and $\mu_b$. Both $V_{ext}$, $\psi_a$ and $\psi_b$ are defined on the domain $x\in[0,L]$. Functions $\psi_a(x)$ and $\psi_b(x)$ are complex functions. The boundary conditions are the periodic ones, i.e.:
$$ \psi_a(x+L)=\psi_a(x) \qquad \psi_b(x+L)=\psi_b(x) $$
Notice that the period of the eigenfunctions should be $L$ while the period of the external potential is $L/3$.
Eventually, there is a constraint on the norms of $\psi_a$ and $\psi_b$, namely: $$ \int_0^L |\psi_a|^2 \, \mathrm{d}x= N, \qquad \int_0^L |\psi_b|^2 \, \mathrm{d}x= M $$
Can you please give me a good strategy to handle this problem and, possibly, a Matlab code?
