Solution (source code)

= Solution

Since each species number $N_i=\int|\psi_i|^2d^3x$ is separately conserved, minimize $E-\mu_1N_1-\mu_2N_2$ with two independent <Lagrange multipliers>. Variation with respect to each complex-conjugate field, with boundary terms removed, gives
$$
\boxed{\left[-\frac{\hbar^2}{2m_1}\nabla^2+V_1+U_{11}|\psi_1|^2+U_{12}|\psi_2|^2\right]\psi_1=\mu_1\psi_1,}
$$
$$
\boxed{\left[-\frac{\hbar^2}{2m_2}\nabla^2+V_2+U_{22}|\psi_2|^2+U_{12}|\psi_1|^2\right]\psi_2=\mu_2\psi_2.}
$$
These are the two stationary <Gross–Pitaevskii equations>; the time-dependent equilibrium fields are $e^{-i\mu_i t/\hbar}\psi_i(\mathbf x)$. The cross term contributes $U_{12}|\psi_j|^2\psi_i$ to each equation, with no factor two, while differentiating the self-interaction term cancels its factor $1/2$. Separate number conservation does not require the two <chemical potentials> to be equal.