Is $D$ the distance of the point charge from the centre of the sphere, or the distance of the point charge from the closest point on the sphere? I'm going to assume the former.
This can be thought of as a special case of a more general problem, where we have
- the conducting sphere of radius $R$, as in your problem, and
- a non-conducting spherical shell of radius $D$, with $\theta$-dependent charge distribution glued into place, which we'll call $\sigma(\theta)$.
In this case, the "charge distribution" has all the charge concentrated at the point $\theta = 0$. We'll discuss later how to model this mathematically.
I will assume that the OP is already familiar with solving Laplace's equation with separation of variables, and has fully worked some problems using this technique in the past. Therefore, I will skip over some of the details.
The first step is always to identify a set of charge-free regions (so that Laplace's equation holds in each region). Since there is a charged shell at $r = D$, the charge-free regions should be
- Region I: $R < r < D$
- Region II: $r > D$
The region $r < R$ is not of interest here; it is known that the spherical shell will completely screen out any electric field from the point charge, so the potential in this region will be the same as the potential on the surface $r = R$.
In Region I and II, respectively, write
\begin{align}
V_I(r, \theta) &= \sum_{\ell=0}^\infty A_\ell r^\ell P_\ell(\cos \theta) + B_\ell r^{-\ell-1} P_\ell(\cos \theta) \\
V_{II}(r, \theta) &= \sum_{\ell=0}^\infty C_\ell r^\ell P_\ell(\cos\theta) + F_\ell r^{-\ell-1} P_\ell(\cos\theta)
\end{align}
The boundary conditions are:
- In Region I, the potential approaches a constant at $r = R$, but we don't know what that constant is yet. That is, at $r = R$, there's no dependence on $\theta$. Since all $P_\ell$'s carry $\theta$ dependence except when $\ell =0$, this leads to the conclusion that $A_\ell r^\ell + B_\ell r^{-\ell-1}$ must vanish whenever $\ell > 0$ and $r = R$.
- In Region II, the potential must go to zero at infinity. This implies that all $C_\ell$ coefficients vanish.
- In Region II, the potential at long distances is dominated by the $1/r$ term, and the coefficient of that term depends only on the total charge of the system, which is $Q$. In particular the potential should go as $V_{II}(r) \approx Q/(4\pi r)$. (I am using units where $\epsilon_0 = 1$ for now, but we'll put it back in later.) That is, $F_0 = Q/(4\pi)$.
- At the boundary between Regions I and II, namely at $r = D$, the potential $V$ will be continuous. (In the specific case of a point charge, it won't be continuous, because it diverges at $(r, \theta) = (D, 0)$, but we can solve the more general case and then take the limit as the charge distribution approaches that of a point charge.) The continuity condition is $V_I(D, \theta) = V_{II}(D, \theta)$ for all $\theta$, and this can hold for all $\theta$ only when the coefficients of $P_\ell(\cos\theta)$ are equal.
- Finally, the areal charge density $\sigma(\theta)$ at $r = D$ leads to a discontinuity in the transverse component of the electric field. In the case of a spherical shell of charge, the transverse component of the electric field is the radial component. The boundary condition derived from Gauss's law is that $E_+ - E_- = \sigma$, where $E_+$ is the radial component of the electric field immediately outside the shell, $E_-$ is the radial component of the electric field immediately inside the shell, and $\sigma$ is the areal charge density at the observation point. Note that $E_+ = -\partial V_{II}/\partial r$ and $E_- = -\partial V_{I}/\partial r$.
I'll skip over most of the algebra, but the final boundary condition leads to:
\begin{align}
(Q/(4\pi) - B_0) D^{-2} + \sum_{\ell=1}^\infty P_\ell(\cos\theta)[(\ell+1) F_\ell D^{-\ell-2} + A_\ell(\ell D^{\ell-1} - (\ell+1) R^{2\ell+1} D^{-\ell-2})] = \sigma(\theta)
\end{align}
So, all we have to do is expand $\sigma(\theta)$ in terms of the spherical harmonics $P_\ell(\cos\theta)$, i.e.
$$ \sigma(\theta) = \sum_{\ell=0}^\infty G_\ell P_\ell(\cos\theta) $$
whereupon we can equate coefficients of $P_\ell(\cos\theta)$ on the left and right sides of the previous equation. Assuming that the $G_\ell$'s are all known, we can combine all the boundary conditions and do a bit of algebra to obtain:
\begin{align}
A_0 &= G_0 D \\
B_0 &= \frac{Q}{4\pi} - G_0 D^2 \\
A_\ell &= \frac{G_\ell}{(2\ell+1) D^{\ell-1}} \qquad (\ell > 0) \\
F_\ell &= \frac{G_\ell}{2\ell+1} (D^{\ell+2} - R^{2\ell+1}/D^{\ell-1}) \qquad (\ell > 0)
\end{align}
Also, note that only $G_0$ contributes to the total charge of the shell; the higher $G_\ell P_\ell(\cos\theta)$ terms integrate to 0. So $G_0 = Q/(4\pi D^2)$. Finally we can write the potentials in regions I and II:
\begin{align}
V_I(r,\theta) &= G_0 D + \sum_{\ell=1}^\infty \frac{G_\ell}{(2\ell+1)D^{\ell-1}}(r^\ell + R^{2\ell+1} r^{-\ell-1}) P_\ell(\cos\theta) \\
V_{II}(r,\theta) &= \frac{G_0 D^2}{r} + \sum_{\ell=1}^\infty \frac{G_\ell}{2\ell+1} (D^{\ell+2} - R^{2\ell+1}/D^{\ell-1}) r^{-\ell-1} P_\ell(\cos\theta)
\end{align}
Now we can solve the specific problem of a point charge, which we'll assume to be located on the positive z-axis. What is the distribution $\sigma(\theta)$ corresponding to a point charge? Well, it's a bit tricky because it must be the case that when it's integrated on the entire surface of the sphere $r = D$, we get $Q$, the total charge. That is
$$ \int_0^\pi \int_0^{2\pi} \sigma(\theta) (D^2 \sin\theta \, \mathrm{d}\varphi \, \mathrm{d}\theta) = Q $$
Even if $\sigma(\theta)$ contains a delta function, the integral will still be zero because $\sin\theta = 0$ when $\theta = 0$. So in fact what we need here is something "stronger" than the delta function, namely its derivative. It turns out that
$$ \sigma(\theta) = -\frac{Q}{2\pi D^2} \delta'(\theta) $$
But the really important question is not what $\sigma$ is, but how to express it as an infinite series in the Legendre polynomials $P_\ell(\cos\theta)$. In other words, what are the coefficients $G$ such that the following holds?
$$ \sum_{\ell=0}^\infty G_\ell P_\ell(\cos\theta) = -\frac{Q}{2\pi D^2} \delta'(\theta) $$
To calculate the $G_\ell$'s, we use the orthogonality relation $\int_0^\pi P_\ell(\cos\theta) P_{\ell'}(\cos\theta) \sin\theta \, \mathrm{d}\theta = \frac{2}{2\ell+1} \delta_{\ell\ell'}$. So
$$ G_\ell = -\frac{2\ell+1}{2} \frac{Q}{2\pi D^2} \int_0^\pi P_\ell(\cos\theta) \sin\theta \, \delta'(\theta) \, \mathrm{d}\theta $$
Using integration by parts, we find that the value of the integral is -1 times the derivative of $P_\ell(\cos\theta)\sin\theta$ evaluated at $\theta=0$. To evaluate that derivative, we can use the product rule. One of the terms has $\sin\theta$ left undifferentiated, evaluating to zero. The other term is $P_\ell(\cos\theta)$ evaluated at $\theta=0$, multiplied by the derivative of $\sin\theta$ at $\theta=0$, which is just $\cos 0 = 1$. And since $P_\ell(\cos\theta)$ at $\theta=0$ is just $P_\ell(1)$, which is 1, we conclude that the entire integral is -1. Therefore
$$ G_\ell = \frac{(2\ell+1)Q}{4\pi D^2} $$
And substituting these values of $G_\ell$ back into the general solution for arbitrary $\sigma$, and putting the factor of $1/\epsilon_0$ back in, we finally obtain
\begin{align}
V_I(r,\theta) &= \frac{Q}{4\pi\epsilon_0} \left[ \frac{1}{D} + \sum_{\ell=1}^\infty \frac{r^\ell + R^{2\ell+1} r^{-\ell-1}}{D^{\ell+1}} P_\ell(\cos\theta) \right] \\
V_{II}(r,\theta) &= \frac{Q}{4\pi\epsilon_0} \left[\frac{1}{r} + \sum_{\ell=1}^\infty (D^\ell - R^{2\ell+1}/D^{\ell+1}) r^{-\ell-1} P_\ell(\cos\theta) \right]
\end{align}
Showing that this is equal to the result obtained from the method of images is left as an exercise to the reader.
By the way, did you notice that each $G_\ell$ coefficient contributes only to the term in the potential that contains $P_\ell$? This is not an accident. It follows from the fact that the differential operator $\nabla^2$ can be separated into a radial part and an angular part, and the angular part, up to scaling, is just $L^2$. And the spherical harmonics are just simultaneous eigenfunctions of $L^2$ and $L_z$, so $\nabla^2$ can only map each spherical harmonic to the same spherical harmonic, with some scaling depending on the radial coordinate and the value of $\ell$.
For example, when a spherical shell's charge distribution is of the form $\cos\theta$, then the potential produced by that shell can only have angular dependence of the form $\cos\theta$ (but obviously still depends on the radial coordinate).
Using this fact, we can derive the general solution above (i.e. for arbitrary $\sigma(\theta)$) much more easily. The idea is that the charge distribution that is induced on the conducting sphere must be such that, for all points in the sphere's interior, every component of the interior multipole moment that is produced by the point charge outside the sphere is completely cancelled by the corresponding interior multipole moment that develops on the surface of the sphere. (The sole exception is, of course, the (0, 0) component that is proportional to the total charge; a neutral conducting sphere cannot screen this out.) This lets us solve for the charge distribution that is induced on the conducting sphere, as an expansion in spherical harmonics, where each coefficient only depends on the corresponding coefficient of the spherical harmonic expansion for $\sigma$. Finally, at each radial coordinate of interest, we can write the potential as an expansion of spherical harmonics, where the coefficient of each term depends only on the coefficients of the same spherical harmonic in $\sigma$ and in the induced charge.
This approach also lets us handle $\varphi$-dependent charge distributions in with only a minimum of additional pain, whereas the explicit separation of variables approach becomes extremely painful when $\varphi$-dependent potentials must be computed.
