Solution (source code)

= Solution

Using $\rho=\varepsilon_0\nabla\cdot E$, $E=-\nabla\phi$, integration by parts, and $\phi=0$ on the boundary gives
$$
U=\frac{\varepsilon_0}{2}\int_V|E|^2\,d^3x.
$$
Put $C=q/(4\pi\varepsilon_0)$. <Gauss's law> gives the radial field
$$
E_r=\begin{cases}0,&r<R,\\C/r^2,&R<r<2R,\\-C/r^2,&2R<r<3R,\\0,&r>3R.\end{cases}
$$
Taking zero potential outside,
$$
\phi=\begin{cases}C/(3R),&r<R,\\C(1/r-2/(3R)),&R<r<2R,\\C(1/(3R)-1/r),&2R<r<3R,\\0,&r>3R.\end{cases}
$$
The charge formula gives
$$
U=\frac12\sum_iQ_i\phi(r_i)=\frac{q^2}{12\pi\varepsilon_0R}.
$$
The field formula gives the same result after integrating $4\pi r^2E_r^2$ over the two nonzero annuli.