Actually, I feel rather retarded for posting this question. Below is the answer.
The equation in question is indeed the Schrödinger's equation for a standard harmonic oscillator. See this link.
It can be written as
$$-\frac{1}{2}\frac{\partial^2 \phi_1(x^1)}{(\partial x^1)^2} + \frac{\omega^2}{2} \left(x^1-\frac{k_2}{eB} \right)^2 \phi_1(x^1) = E' \phi_1(x^1),$$
where $\omega = e B$ and $E' = \frac{1}{2}(E^2 - (k_3-eA)^2 + eB)$.
$\frac{k_2}{eB}$ simply shifts the center of oscillation.
We have $E' = \left(\frac{1}{2} + n_1 \right)\omega$. Therefore
$$E^2 = (k_3-eA)^2 + 2 n_1 eB = \left(\frac{2\pi n_3}{L}-eA\right)^2 + 2 n_1 eB.$$
This actually looks better than (19.9) in Zhong-Zhi Xianyu's solution.
However, I have not figured out part (e).
As far as the equation above is concerned, when a certain $n_3$ become invalid as we change background $A$ adiabatically, $-n_3$ at the same time would become valid. ($n_3$ can be either positive or negative integers.)
My guess is that the first-order equations of $\phi_1$ and $\phi_2$ will put further restrictions on what values of $n_3$ are actually allowed.