To be honest, I just learned about all this myself in the last months, so I am not sure whether this is actually correct.
Since you have this spherical symmetry, I think that you need spherical harmonics. They are orthogonal functions, think of them as a Fourier series on the surface of a sphere.
Your charge density $\sigma$ does not depend on $\phi$, therefore we can use the simpler Legendre-polynomials, where $m = 0$. Fist, we want to express the potential like so:
$$
\varphi(\vec x) = \frac1{\varepsilon_0} \sum_{l=0}^\infty \frac{1}{2l+1} \varrho_{l} \frac{1}{r^{l+1}} Y_{l,0} (\theta, \phi)
$$
The coefficients $\varrho_l$ are given by:
$$
\varrho_l = \int \mathrm d^3 x'\, r'^l \varrho(\vec x) Y^*_{l,0} (\theta', \phi')
$$
Now with
$$
Y^*_{l,0}(\theta', \phi') = \sqrt{\frac{2l+1}{4\pi}} P_l(\cos \theta')
$$
and
$$
P_0(x) = 1
,\quad
P_1(x) = x
$$
we can calculate the coefficitents.
But first, we need to convert the surface charge density $\sigma$ into a volume charge density. For that, we use the $\delta$-distribution:
$$
\varrho(\vec x) = \delta(r-R) \sigma_0 \cos(\theta)
$$
If you plug those into the $\varrho_l$, you will get $\varrho_0 = 0$ and $\varrho_1 = \sqrt{\frac{3}{4 \pi}} \frac{4}{3} \pi R^3 \sigma_0$. I hope this is correct.
Then we can put this into the first formula and get $\varphi$:
$$
\varphi(\vec x) = \frac{1}{\varepsilon_0} \sqrt{\frac{3}{4 \pi}} \frac{4}{3} \pi R^3 \frac{1}{r^2} \sigma_0 \cos \theta
$$
Since this is a pure dipol potential, the $1/r^2$ seems about right. And if you look at the dimensions, the $R^3 \sigma_0 / r^2$ have just the needed Charge/Length.
Spherical harmonics might be overkill, maybe there is a simpler method to do this.