Suppose a concentric cylinder of height $c$ and radius $b$ with it's top cap held at $V=V_0$ and all other surfaces held at $V=0$. For calculating the potential inside the cylinder I use the Laplace equation because there is no free charge. For solving the Laplace equation in cylindrical coordinates I use a product approach:
\begin{align*}
\varphi(r,z,\phi)=\mathcal{R}(r) \cdot \mathcal{Z}(z) \cdot \Psi(\phi)
\end{align*}
With the followings solutions:
\begin{align*}
\mathcal{Z} &\propto e^{\pm \, k z} \\
\Psi &\propto e^{\pm \, i \nu \phi} \\
\mathcal{R} &\propto J_{\nu}(k r) \ \mathrm{or} \ N_{\nu}(k r)
\end{align*}
Where $J_{\nu}$ are the Bessel functions and $N_{\nu}$ are the Neumann functions. Our boundary conditions are:
\begin{align*}
\varphi(r,z=c)&=V_0 \\
\varphi(r,z=0)&=0 \\
\varphi(r=b,z)&=0
\end{align*}
Because of rotational symmetry we can find $\nu=0$ and therefore we can neglect the $\Psi$ component. Now following the calculation from the following lecture on page 4 we end up with the solution:
\begin{align*}
\varphi(r,z) = \sum_{n=0}^{\infty} A_n \, J_0(k_n r) \, \sinh(k_n z) \\
A_n = \frac{V_0}{\frac{b^2}{2} [J_1(k_n b)]^2 \, \sinh(k_n c)} \int_0^b r \, J_0(k_n r) \, dr
\end{align*}
($A_n$ is obtained by making use of the orthogonality of the Bessel functions, see page 3.) I have already checked this solution numerically and it's correct.
So far so good. But now suppose the same cylinder, filled with a dielectric material of $\epsilon=\epsilon_1$ from $r=0$ to $r=a$ and another dielectric material with $\epsilon=\epsilon_2$ from $r=a$ to $r=b$. Now the problem becomes more complicated and we get an additional boundary condition:
\begin{align*}
\varphi(r,z=c)&=V_0 \quad \quad (1) \\
\varphi(r,z=0)&=0 \quad \quad (2) \\
\varphi(r=b,z)&=0 \quad \quad (3) \\
\varphi(r=a,z)&=f(z) \quad \quad (4) \\
\end{align*}
To solve this problem I have tried two different approaches with yet no success. In both approaches I divided the solution into two solution. $\varphi_1$ in the area of $\epsilon_1$ and $\varphi_2$ in the area of $\epsilon_2$. First approach goes as follows:
Approach 1
$\Psi_{1,2}$: Rotational symmetry $\rightarrow$ $\nu=0$ and $\Psi=1$
$\mathcal{Z}_{1,2}$: Due to boundary condition (1) and (2) $\mathcal{Z}=\sinh(k_{1n} z)$
$\mathcal{R}_{1 \ }$: Because $N_0(0) = \infty$ , $\mathcal{R}_{1}=J_0(k_{1n} r)$
$\mathcal{R}_{2 \ }$: We combine Bessel and Neumann functions so that that they hold boundary condition (3) $\rightarrow$ $\mathcal{R}_{2} = G_0 = \frac{J_0(k_{1n} r)}{J_0(k_{1n} b)} - \frac{N_0(k_{1n} r)}{N_0(k_{1n} b)}$
The resulting solutions are:
\begin{align*}
\varphi_1 &= \sum_{n=0}^{\infty} A_{1n} \, \sinh(k_{1n} z) \, J_0(k_{1n} r) \\
\varphi_2 &= \sum_{n=0}^{\infty} A_{2n} \, \sinh(k_{2n} z) \, G_0(k_{2n} r) \\
\end{align*}
with
\begin{align*}
A_{1n} &= \frac{V_0}{\sinh(k_{1n} c) \frac{a^2}{2} \left( [J_0(k_{1n} c)]^2 + [J_1(k_{1n} c)]^2 \right)} \int_0^a r \, J_0(k_n r) \, dr \\
A_{2n} &= \frac{1}{\sinh(k_{2n} c) \int_a^b r \, G_0(k_{2n} r) \, dr} \int_a^b r \, G_0(k_{2n} r) \, V_0 \, dr = \frac{V_0}{\sinh(k_{2n} c)}
\end{align*}
Now the only two parameters undetermined are $k_{1n}$ and $k_{2n}$. To solve for the two unknowns we have the two continuity conditions:
\begin{align*}
\varphi_1(r=a,z) &= \varphi_2(r=a,z) = f(z) \\
\epsilon_1 \frac{\partial \varphi_1(r,z)}{\partial r} \biggr\rvert_{r=a} &= \epsilon_2 \frac{\partial \varphi_2(r,z)}{\partial r} \biggr\rvert_{r=a}
\end{align*}
But it turns out, that the equations are just too complicated to solve for $k_{1n}$ and $k_{2n}$. So I got stuck at this point here.
Approach 2
So for me the more promising approach is to use the superposition principle and devide the problem in one problem A with the boundary conditions:
\begin{align*}
\varphi(r,z=c)&=V_0 \quad \ (1A) \\
\varphi(r,z=0)&=0 \quad \quad (2A) \\
\varphi(r=b,z)&=0 \quad \quad (3A) \\
\varphi(r=a,z)&=0 \quad \quad (4A) \\
\end{align*}
and another problem B with the boundary conditions:
\begin{align*}
\varphi(r,z=c)&=0 \quad \quad (1B) \\
\varphi(r,z=0)&=0 \quad \quad (2B) \\
\varphi(r=b,z)&=0 \quad \quad (3B) \\
\varphi(r=a,z)&=f(z) \ \ (4B) \\
\end{align*}
If I solve both of these problems and add up their solutions, due to the superposition principle I get the solution I desire:
\begin{align*}
\varphi = \varphi_A + \varphi_B
\end{align*}
The solution to $\varphi_A$ is already given in the lecture I have posted above, see link on page 6-7. For $\varphi_B$ it goes as follows:
$\Psi^B_{1,2}$: Rotational symmetry $\rightarrow$ $\nu=0$ and $\Psi=1$
$\mathcal{Z}^B_{1,2}$: Because neither $\sinh$ nor $\cosh$ can satisfy $\varphi(r,z=c)=0$ and $\varphi(r,z=0)=0$, $k$ has to become complex $k \rightarrow ik$ so that the $\sinh$ will become a $\sin$ with $k_n= i \frac{n \pi}{c}$: $\mathcal{Z}^B_{1,2}= \sin(\frac{n \pi}{c} z)$
$\mathcal{R}^B_{1 \ }$: Similar to the first approch we get: $\mathcal{R}^B_{1} = J_0(i \frac{n \pi}{c} r)$
$\mathcal{R}^B_{2 \ }$: And $\mathcal{R}^B_{1} = G_0 = \frac{J_0(i \frac{n \pi}{c} r)}{J_0(i \frac{n \pi}{c} b)} - \frac{N_0(i \frac{n \pi}{c} r)}{N_0(i \frac{n \pi}{c} b)}$
The resulting solutions are:
\begin{align*}
\varphi_A = \sum_{n=0}^{\infty} A_{n} \, \sinh ( k_n z) \, g_0 (k_n r)
\end{align*}
with
\begin{align*}
A_n = \frac{V_0}{\sinh(k_n c)} \quad \quad \mathrm{and} \quad \quad g_0 = \frac{J_0(k_n r)}{J_0(k_n a)} - \frac{N_0(k_n r)}{N_0(k_n a)}
\end{align*}
and
\begin{align*}
\varphi_{B,1} &= \sum_{n=0}^{\infty} B_{1n} \, \sin \left( \frac{n \pi}{b} z \right) \, J_0 \left(i \frac{n \pi}{c} r \right) \\
\varphi_{B,2} &= \sum_{n=0}^{\infty} B_{2n} \, \sin \left( \frac{n \pi}{b} z \right) \, G_0 \left(i \frac{n \pi}{c} r \right)
\end{align*}
with
\begin{align*}
B_{1n} &= \frac{1}{\frac{c}{2} J_0 \left(i \frac{n \pi}{c} a \right)} \int_0^c f(z) \sin \left(\frac{n \pi}{b} z \right) \, dz \\
B_{2n} &= \frac{1}{\frac{c}{2} G_0 \left(i \frac{n \pi}{c} a \right)} \int_0^c f(z) \sin \left(\frac{n \pi}{b} z \right) \, dz
\end{align*}
The big advantage to this approach is that $k_n$ is easily determined by finding the $k_n$ for which $g_0(b)=0$. But the big disadvantage is that there is no way, to get the $f(z)$ out of the definite integral in $B_{1n}$ and $B_{2n}$. Here $f(z)$ is the only unknown left. But by using the continuity condition I see no way to find out $f(x)$. From my numerical solution and also in approach 1, one can see, that $f(x)$ has to be a $sinh(...)$. Maybe this is a hint for you.
I think there is no way to solve the equation by using the definite integral of $f(x)$ to express $B_{1n}$ and $B_{2n}$.
In Germany we have a saying "...den Wald vor lauter Bäumen nicht sehen." -- "...overlooking the forest due to too many trees.". And I think, and I hope that it is something obvious that I'm overlooking and I hope you can help me to find it. Also I am open for new proposals to solve this equation.
Big thanks in advance!
Pearli