Let be
$$\frac{2a}{Q}V(\theta,\varphi)=f(\theta,\varphi)=2\sin\theta\cos\varphi+\cos^2\theta.\tag 1$$
The Laplace spherical harmonics form a complete set of orthonormal functions and thus form an orthonormal basis of the Hilbert space of square-integrable functions. On the unit sphere, any square-integrable function can thus be expanded as a linear combination of these:
$$
f(\theta,\varphi)=\sum_{\ell=0}^\infty \sum_{m=-\ell}^\ell f_\ell^m \, Y_\ell^m(\theta,\varphi)\tag 2
$$
where $Y_\ell^m( \theta , \varphi )$ are the Laplace spherical harmonics defined as
$$
Y_\ell^m( \theta , \varphi ) = \sqrt{{(2\ell+1)\over 4\pi}{(\ell-m)!\over (\ell+m)!}} \, P_\ell^m ( \cos{\theta} ) \operatorname{e}^{i m \varphi } =N_{\ell}^m P_\ell^m ( \cos{\theta} ) \operatorname{e}^{i m \varphi }\tag 3
$$
and where $N_{\ell}^m$ denotes the normalization constant
$ N_{\ell}^m \equiv \sqrt{{(2\ell+1)\over 4\pi}{(\ell-m)!\over (\ell+m)!}},$ and $P_\ell^n(\cos\theta)$ are the associated Legendre polynomials.
The Laplace spherical harmonics are orthonormal
$$
\int_{\theta=0}^\pi\int_{\varphi=0}^{2\pi}Y_\ell^m \, Y_{\ell'}^{m'*} \, d\Omega=\delta_{\ell\ell'}\, \delta_{mm'},
$$
where $δ_{ij}$ is the Kronecker delta and $\operatorname{d}\Omega = \sin\theta \operatorname{d}\varphi\operatorname{d}\theta$.
The expansion coefficients are the analogs of Fourier coefficients, and can be obtained by multiplying the above equation by the complex conjugate of a spherical harmonic, integrating over the solid angle $Ω$, and utilizing the orthogonality relationships. This is justified rigorously by basic Hilbert space theory. For the case of orthonormalized harmonics, this gives:
$$
f_\ell^m=\int_{\Omega} f(\theta,\varphi)\, Y_\ell^{m*}(\theta,\varphi)\operatorname{d}\Omega = \int_0^{2\pi}\operatorname{d}\varphi\int_0^\pi \operatorname{d}\theta\,\sin\theta f(\theta,\varphi)Y_\ell^{m*} (\theta,\varphi). \tag 4
$$
where $ Y_\ell^{m*} (\theta, \varphi) = (-1)^m Y_\ell^{-m} (\theta, \varphi)$.
The evaluation of the expansion $f_\ell^m$ may be very long in this way...
We can use some tricks in your case for $f(\theta,\varphi)=2\sin\theta\cos\varphi+\cos^2\theta$ observing that
$$
\sin\theta=-P_1^1(\cos\theta)
$$
and
$$
Y_{1}^{-1}(\theta,\varphi) - Y_{1}^{1}(\theta,\varphi)= {1\over 2}\sqrt{3\over 2\pi}\cdot e^{-i\varphi}\cdot\sin\theta-{-1\over 2}\sqrt{3\over 2\pi}\cdot e^{i\varphi}\cdot\sin\theta =\sqrt{3\over 2\pi} \sin\theta\cos\phi
$$
so that
$$
2\sin\theta\cos\varphi=2\sqrt{\frac{2\pi}{3}}\left(Y_1^{-1}(\theta,\varphi)-Y_1^{1}(\theta,\varphi)\right)=-2P_1^1(\cos\theta)\cos\varphi\tag 5
$$
and observing that
$$
\cos^2\theta=\frac{1}{3}P_0^0+\frac{2}{3}P_2^0
$$
and using the relation $Y_\ell^0(\theta,\varphi)=\sqrt{\frac{2\ell+1}{4\pi}}P_\ell^0(\cos\theta)$ where $P_\ell^0(\cos\theta)$ are the ordinary Legendre's polynomials $P_\ell(\cos\theta)$, we obtain
$$
\cos^2\theta=\frac{1}{3}P_0^0(\cos\theta)+\frac{2}{3}P_2^0(\cos\theta)=2\sqrt{\pi}Y_0^0(\theta,\varphi)+\frac{4}{3}\sqrt{\frac{\pi}{5}}Y_2^0(\theta,\varphi).\tag 6
$$
Finally, putting together (5) and (6) in (1) we obtain
$$
f(\theta,\varphi)=2\sqrt{\frac{2\pi}{3}}Y_1^{-1}(\theta,\varphi)-2\sqrt{\frac{2\pi}{3}}Y_1^{1}(\theta,\varphi)+2\sqrt{\pi}Y_0^0(\theta,\varphi)+\frac{4}{3}\sqrt{\frac{\pi}{5}}Y_2^0(\theta,\varphi)\tag 7
$$
so that, comparing (7) and (2), the coefficients $f_\ell^m$ are
$$
f_1^{-1}=2\sqrt{\frac{2\pi}{3}}\qquad
f_1^{1}=-2\sqrt{\frac{2\pi}{3}}\qquad
f_0^{0}=2\sqrt{\pi}\qquad
f_2^{0}=\frac{4}{3}\sqrt{\frac{\pi}{5}}\tag 8
$$
So you have
$$\small
V(\theta,\varphi)=\frac{Q}{a}\left[\sqrt{\frac{2\pi}{3}}Y_1^{-1}(\theta,\varphi)-\sqrt{\frac{2\pi}{3}}Y_1^{1}(\theta,\varphi)+\sqrt{\pi}Y_0^0(\theta,\varphi)+\frac{2}{3}\sqrt{\frac{\pi}{5}}Y_2^0(\theta,\varphi)\right]
$$
The general solution to the Laplace equation outside the sphere is
$$
\Phi(r,\theta,\varphi)=\sum_{\ell=0}^\infty \sum_{m=-\ell}^\ell \frac{B_\ell^m}{r^{\ell+1}} \, Y_\ell^m(\theta,\varphi)\tag 9
$$
with $B_\ell^m=\frac{Q}{a}f_\ell^m$,
that is
$$\small
\Phi(r,\theta,\varphi)=\frac{Q}{a}\left[\frac{1}{r^2}\left(\sqrt{\frac{2\pi}{3}}Y_1^{-1}(\theta,\varphi)-\sqrt{\frac{2\pi}{3}}Y_1^{1}(\theta,\varphi)\right)+\frac{\sqrt{\pi}}{r}Y_0^0(\theta,\varphi)+\frac{1}{r^3}\left(\frac{2}{3}\sqrt{\frac{\pi}{5}}\right)Y_2^0(\theta,\varphi)\right]
$$
Best Answer
The problem is radially symmetric both in and outside the sphere. Because of the radial symmetry only the component with $l=0$ and $m=0$ survives; so that the potential is only a function of $r$.
As you correctly noted, the potential thus satisfies the equation $$ \frac{1}{r^2}\frac{d}{d r}\left(r^2 \frac{d V(r)}{d r}\right)=-\rho(r)$$ with $$\rho(r) = \begin{cases} \rho_0 r, & r<R, \\ 0, &R>0.\end{cases}$$
Solving the potential outside the sphere $$\frac{1}{r^2}\frac{d}{d r}\left(r^2 \frac{d V(r)}{d r}\right)= 0,$$ we find $V(r>R) = \frac{C}{r}$ with $C$ a constant; we set the integration constant such that $V$ approaches $0$ for $r\to\infty$.
The equation for the potential inside the sphere fulfils the equation $$\frac{1}{r^2}\frac{d}{d r}\left(r^2 \frac{d V(r)}{d r}\right)= -\rho_0 r$$ with the solution $$V(r<R) = -\frac{\rho_0 r^3}{12} + c $$ where we took into account that the potential should remain finite for $r=0$.
From the condition that $V(r)$ and $V'(r)$ are continuous at $r=R$, we obtain the relations $$ c- \frac{\rho_0 R^3}{12} = \frac{C}R, \qquad -\frac{\rho_0 R^2}4 = -\frac{C}{R^2}$$ with the solution $C= \rho_0 R^4/4$ and $c= \rho_0 R^3/3$.
Thus the solution has the form $$ V(r) = \begin{cases} \frac{\rho_0}{12} (4 R^3 -r^3),& r<R,\\ \frac{\rho_0 R^4}r, &r<R. \end{cases}$$