I'd say the most intuitive way is to integrate the electrostatic equations (written in Heaviside-Lorentz units and for the vacuum),
[tex]\vec{\nabla} \times \vec{E}=0, \quad \vec{\nabla} \cdot \vec{E}=\rho.[/tex]
The first equation tells us that the electric field has a scalar potential,
[tex]\vec{E}=-\vec{\nabla} \cdot \Phi.[/tex]
Gauß's Law then becomes the Poisson Equation,
[tex]\Delta \Phi=-\rho.[/tex]
In your case of a homogeneously charged spherical shell you have
[tex]\rho(\vec{x})=\frac{Q}{4 \pi R^2} \delta(r-R),[/tex]
where [itex]r=|\vec{x}|[/itex].
The charge distribution is spherically symmetric and thus one can expect that the potential is a function of [itex]r[/itex] only. Writing the Laplace operator in terms of spherical coordinates our Poisson equation reads
[tex]\frac{1}{r} \frac{\mathrm{d}^2}{\mathrm{d} r^2} (r \Phi)=-\frac{Q}{4 \pi R^2} \delta(r-R). \qquad (*)[/tex]
For [itex]r \neq R[/itex] the right-hand side is 0 and the general solution thus reads
with integration constants [itex]A[/itex] and [itex]B[/itex], which are to be chosen different for [itex]r<R[/itex] and [itex]r>R[/itex] due to the singularity of the [itex]\delta[/itex] distribution. We have to find the constants from appropriate boundary conditions.
Start with [itex]r<R[/itex]. Since there are no sources or other singularities there, we must have [itex]A=0[/itex] in that region, because otherwise the potential would have a singularity in the origin. Thus we have
[tex]\Phi(r)=B \quad \text{for} \quad r<R.[/tex]
For [itex]r>R[/itex] we only demand that the potential vanishes at infinity. This fixes simply the overall additive constant of the potential, which has no physical meaning anyway.
So our solution reads
[tex]\Phi(r)=B \Theta(R-r)+\frac{A}{r} \Theta(r-R).[/tex]
The potential should be continuous at [itex]r=R[/itex] since only its second derivative should have a [itex]\delta[/itex]-distribution like singularity. The first derivative must thus have a jump. So we must have [itex]B=A/R[/itex]. And thus
[tex]\Phi(r)=A \left [\frac{1}{R} \Theta(R-r)+\frac{1}{r} \Theta(r-R) \right ].[/tex]
To determine finally [itex]A[/itex], we have to multiply the equation (*) by [itex]r[/itex] and integrate over an infinitesimal interval around [itex]r=R[/itex]. This gives
[tex]\left ([r \Phi(r)]' \right )_{r=R+0^+}-\left ([r \Phi(r)]' \right )_{r=R-0^+}=-\frac{Q}{4 \pi R}.[/tex]
Evaluating the left-hand side gives
[tex]-A/R=-\frac{Q}{4 \pi R} \; \Rightarrow \; A=\frac{Q}{4 \pi}.[/tex]
The potential thus reads
[tex]\Phi(r)=\frac{Q}{4 \pi} \left [\frac{\Theta(R-r)}{R}+\frac{\Theta(r-R)}{r} \right ].[/tex]
Inside the spherical shell the potential is constant and thus the field vanishes there. Outside you have a Coulomb potential as if the charge is sitting at the center.
Of course you can also evaluate the potential using the Green's-function method, according to which you just integrate over all points of your charge distribution, using Coulomb's Law, i.e.,
[tex]\Phi(\vec{x})=\frac{1}{4 \pi} \int_{V} \mathrm{d}^3 \vec{x}' \frac{\rho(\vec{x}')}{|\vec{x}-\vec{x}'|}.[/tex]
In your case the integral simplifies to an integral along the spherical surface,
[tex]\Phi(\vec{x})=\frac{1}{4 \pi} \frac{Q}{4 \pi R^2} R^2 \int_0^{2 \pi} \mathrm{d} \varphi' \int_0^{\pi} \mathrm{d} \vartheta' \frac{\sin \vartheta'}{|R \vec{e}_r(\vartheta',\varphi')-\vec{x}|}.[/tex]
Here I've used spherical coordinates [itex](\vartheta',\varphi')[/itex] to parametrize the spherical shell and set the radial unit vector to [itex]\vec{e}_r(\vartheta',var\phi')=(\cos \varphi' \sin \vartheta',\cos \varphi' \sin \vartheta',\cos \vartheta')[/itex]. Now, you can choose the spherical coordinates such that [itex]\vec{x}=r \vec{e}_z[/itex], i.e., that this vector points along the polar axis of the spherical coordinates. Then you have
[tex]|R\vec{e}_r-\vec{x}|^2=(R^2 + r^2 - 2 R r \cos \vartheta').[/tex]
The integral over [itex]\varphi'[/itex] thus just gives a factor [itex]2 \pi[/itex], and you get
[tex]\Phi(\vec{x})=\Phi(r)=\frac{Q}{8 \pi} \int_0^{\pi} \mathrm{d} \vartheta' \frac{\sin \vartheta'}{\sqrt{R^2+r^2-2 R r \cos \vartheta'}}.[/tex]
Now you substitute [itex]u=\cos \vartheta', \quad \mathrm{d} u=-\mathrm{d} \vartheta' \sin \vartheta'[/itex], which gives
[tex]\Phi(r)=\frac{Q}{8 \pi} \int_{-1}^{1} \mathrm{d} \vartheta' \frac{1}{\sqrt{R^2+r^2-2 R r u}}.[/tex]
This you can integrate in closed form, giving
[tex]\Phi(r)=\frac{Q}{8 \pi} \frac{R+r-|R-r|}{r R} = \frac{Q}{4 \pi} \left [\frac{\Theta(R-r)}{R}+ \frac{\Theta(r-R)}{r} \right].[/tex]
This is of course the same result as with the direct calculation of the Poisson equation.