This is not a direct answer to the question asked, but rather a different approach to the problem. It's much easier to derive the monopole harmonics from a geometric consideration of the Hopf map ${\rm Hopf}: S^3\to S^2$, which is the same thing as the Hopf bundle $\pi: {\rm SU}(2) \to {\rm SU}(2)/{\rm U}(1)$. Then one immediately sees monopole spherical harmonics coincide with the ${\rm SU}(2)$ represenationrepresentation matrices. Recall that the ordinary spherical harmonics are given in terms of these matrices as $$ Y^L_m(\theta, \phi)= \sqrt{\frac{2L+1}{4\pi}}[D^L_{m,0}(\theta, \phi,\psi)]^* $$ where, on the LHS, $\theta$ $\phi$ are the spherical polar angles on $S^2$. On the RHS $\theta$, $\phi$, $\psi$ are the Euler-angle co-ordinates on $S^3$. The Hopf map simply takes $(\theta, \phi, \psi)$ to $(\theta, \phi)$. Because $$ D^L_{m,n}(\theta, \phi,\psi)= <L,m|\exp(-i\phi J_z)\exp(-i\theta J)\exp(-i \psi J_z)|L,n> $$ the RHS does not depend on $\psi$ when $n=0$.
For a monople of strength $\int B d(Area)= 4\pi \Lambda $, where $4\pi \Lambda $ must be an integer multiple of $2\pi$, the monopole harmonics are
$$
{\mathcal Y}^J_{m, \Lambda} (\theta, \phi,\psi)= \sqrt{\frac{2J+1}{4\pi}}[D^J_{m\Lambda }(\theta,\phi,\psi)]^*
$$
Now both sides depend on $\psi$, so one must choose a $\psi$ for each pair $(\theta, \phi)$. This is just a choice of gauge.
For a derivation and explanantion of all this see page 278 in my online lecture notes https://courses.physics.illinois.edu/phys509/sp2017/bmaster.pdf.