# Numerical solution of two coupled nonlinear eigenvalue problems

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_a(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?