Spherical coordinates are used here; let $r$ be radial coordinate,
$\theta$ the polar angle and $\phi$ the azimuthal angle.
Let $V$ denote the electric scalar potential; the unit of $V$ is volte.
Although the coordinate dependence of $V$ is $V(r,\theta,\phi)$,
let us restrict ourselves to the case where $V$ only depends on $r$ for a while.
Assume that the electric scalar potential is given as
\begin{equation*}
V(r)=k_1\frac{1}{r}+k_2\tag{1}
\end{equation*}
where $k_1$ and $k_2$ are constants, given by,
\begin{align*}
k_1=-\frac{R_2R_1(V_2-V_1)}{R_2-R_1} \\
k_2=\frac{R_2V_2-R_1V_1}{R_2-R_1} \tag{2}
\end{align*}
Equations (1) and (2) satisfy the boundary conditions upon metal surfaces:
\begin{align*}
V(R_1)=V_1 \\
V(R_2)=V_2.\tag{3}
\end{align*}
The electric field in this case has only r-component and given as
\begin{equation*}
E_r(r)=-\mathrm{grad}V(r)=-\frac{k_1}{r^2}\tag{4}
\end{equation*}
Up to this point, we have assumed that the physical quantity depends only on $r$.
From this point on, we will change the equation to depend on other coordinate variables as well.
Electric flux density ($\mathbf{D}$), has r-component:
\begin{align*}
D_r(r,\theta,\phi)&=\epsilon(\theta,\phi) E_r(r)=\epsilon(\theta,\phi)\left(-\frac{\partial V(r)}{\partial r}\right)=-\epsilon(\theta,\phi)\frac{k_1}{r^2} \\
&=\epsilon(\theta,\phi)\frac{R_2R_1(V_2-V_1)}{R_2-R_1}\frac{1}{r^2}\tag{5}
\end{align*}
As a special case, let us focus on the case $r = R_1$. Electric flux density becomes
\begin{align*}
D_r(r=R_1,\theta,\phi) =\epsilon(\theta,\phi)\frac{R_2(V_2-V_1)}{R_1(R_2-R_1)}.\tag{6}
\end{align*}
Total electric charge on $r=R_1$ metal surface is
\begin{align*}
Q=\int_{S}D_r(r=R_1,\theta,\phi)dS
=\int_{S}\left\{\epsilon(\theta,\phi)\frac{R_2(V_2-V_1)}{R_1(R_2-R_1)}\right\}dS.\tag{7}
\end{align*}
where $dS$ is the small area at $r=R_1$ sphere (in unit of $m^2$).
If in the particular case that permitivity is uniform,
i.e.,$\epsilon(\theta,\phi)=\epsilon_1$, then
\begin{align*}
Q=\epsilon_i\frac{R_2(V_2-V_1)}{R_1(R_2-R_1)}4\pi R_1^2
=\epsilon_i\frac{4\pi R_1R_2}{R_2-R_1}(V_2-V_1)\tag{8}
\end{align*}
From this equation, the capacity matrix is found to be:
\begin{equation*}
C=\epsilon_i\frac{4\pi R_1R_2}{R_2-R_1}.\tag{9}
\end{equation*}
For another case that the questioner is focusing on,
the capacity is the following equation,
\begin{equation*}
C=\frac{\epsilon_1+\epsilon_2+\epsilon_3+\epsilon_4}{4}\frac{4\pi R_1R_2}{R_2-R_1}.\tag{10}
\end{equation*}
Now let me justify that equations (1), (2), (4), (5) and (6) are the solution
for the case where the permitivity depends on $\theta,\phi$:
\begin{equation*}
\epsilon=\epsilon(\theta,\phi)\;\;\;\text{if }R_1<r<R_2\tag{11}.
\end{equation*}
The general starting electro-static field equation is as follows.
\begin{align*}
\mathrm{div}\mathbf{D}=\rho,\tag{12} \\
\mathbf{D}=\epsilon\mathbf{E},\tag{13} \\
\mathbf{E}=-\mathrm{grad}V.\tag{14}
\end{align*}
The symbols in (12) to (14) are according to the standard conventions.
The electric free charge density $\rho$ is $0$ in the insulating materials
and $\rho \neq 0$ at the metal surfaces.
Thus within insulationg region, (12) becomes
\begin{align*}
\mathrm{div}\left(\epsilon\mathbf{E}\right)
&=\mathrm{grad}\epsilon\cdot\mathbf{E}+\epsilon\mathrm{div}\mathbf{E} \\
&=\frac{1}{r}\frac{\partial\epsilon}{\partial\theta}E_{\theta}
+\frac{1}{r\sin\theta}\frac{\partial\epsilon}{\partial\phi}E_{\phi}
+\epsilon\left\{
\frac{1}{r^2}\frac{\partial(r^2E_r)}{\partial r}
+\frac{1}{r\sin\theta}\frac{\partial}{\partial\theta}(E_{\theta}\sin\theta)
+\frac{1}{r\sin\theta}\frac{\partial E_{\phi}}{\partial\phi}\right\} \\
&=\frac{1}{r}\frac{\partial\epsilon}{\partial\theta}\left(-\frac{1}{r}\frac{\partial V}{\partial\theta}\right)
+\frac{1}{r\sin\theta}\frac{\partial\epsilon}{\partial\phi}\left(-\frac{1}{r\sin\theta}\frac{\partial V}{\partial\phi}\right) \\
&+\epsilon\left\{
\frac{1}{r^2}\frac{\partial\left(r^2\left(-\frac{\partial V}{\partial r}\right)\right)}{\partial r}
+\frac{1}{r\sin\theta}\frac{\partial}{\partial\theta}\left(\left(-\frac{1}{r}\frac{\partial V}{\partial\theta}\right)\sin\theta\right)
+\frac{1}{r\sin\theta}\frac{\partial}{\partial\phi}\left(-\frac{1}{r\sin\theta}\frac{\partial V}{\partial\phi}\right)\right\}\\
&=0\tag{15}
\end{align*}
Equation (15) is the equation to be solved, although it is complex, with V as the unknown.
Starting from equation (15), trying to find V is difficult,
but substituting (1) into (15) shows that it is one of the solutions.
On the other hand, someone has proved the uniqueness of the solution,
so we can conclude that (1) is the solution of (15).
