Electric Potential Field of Parallel Electrodes within Grounded Shell Any help would be greatly appreciated
I think I need to solve this using boundary potential conditions and Laplace equations.
Here is the system I'm attempting to model:

The red rods are metal electrodes at the same potential ($V_r$) and on the same $x$ plane (offset equally from the $z$ axis). The orange cylinder is at $V=0$. I need the potential field map in 3D Cartesian coordinates, with variables of electrode voltage, electrode radius, electrode offset from the center axis, shell radius, overall system length.
Update 1:
I have simulated the system using FEMM, but in order to incorporate the answer with an applied magnetic field and particle motion, I need the solution in equation form.

Update 2:
I have this reference for a single coaxial electrode (http://www.ittc.ku.edu/~jstiles/220/handouts/Example%20The%20Electorostatic%20Fields%20of%20a%20Coaxial%20Line.pdf) in polar coordinates.
 A: Brief Summary
Numerically, the mean value property of harmonic functions allows you to get an approximate solution to boundary value problems relatively quickly. Often you can improve convergence to a solution with a good initial guess, however, so analytical approaches can still be useful.
Consider the limit of an infinitely long cylinder. There is a simple closed form expression for the potential in the limit where the electrode radius vanishes, and for non-vanishing radius it is possible to use an iterative approach to find a sequence of analytic approximations. 
Detailed Summary
Firstly, as mentioned in the comments, for any harmonic function $\phi(x,y)$ there is a holomorphic function $f(x+iy)$ such that $\phi(x,y)$ is the real-part of $f$. Another crucial property of solutions to Laplace's equation is invariance under conformal mappings. If $\nabla^2\phi(x,y)=0$ and $f$ is a conformal map (loosely), then $\nabla^2\phi(f(x+iy))=0$ as well. This property can be used to simplify the problem somewhat. 
Let $R$ be the radius of the large cylinder, and let $\sigma(z)=Ai\frac{i-z/R}{i+z/R}$ (where $A$ is a real constant). You can check that this map sends the interior of the large cylinder (viewed as a disk in the complex plane $\mathbb C$) to the upper half plane $\mathbb H$. You can also check that it maps circles to other circles, so the electrode cylinders are mapped into two rescaled circles in the upper half plane. Moreover, the centers of the electrodes are mapped to a line parallel to the real axis, and both circles in $\mathbb H$ have equal radius. As promised, the picture:

In the new coordinates, the potential $\phi$ must vanish on the real line $z=\bar z$, and have some constant value $V$ on two circles in the upper half plane. To solve this, we can invoke another standard trick of boundary value problems: the method of image charges. If we extend $\phi$ by simultaneously reflecting it across the real-line while reversing its sign, we get a larger (but more symmetric) boundary value problem.
Now I use the fact that a harmonic function is the real part of a holomorphic function. Let $\psi(x,y)$ be the potential produced by charges on the upper right cylinder (centered at point $z_0$ with respect to the origin, which is set at the 'center of gravity' of the augmented problem), and let $f(z-z_0)$ be the holomorphic function that $\psi$ is the real part of ($f$ is viewed as a function defined near the origin but displaced by $z_0$). Then "by symmetry", the fields generated by the other cylinders are (proceeding counter-clockwise) $\overline{f(-\bar z-z_0)}=\bar f(-z-\bar z_0)$, $-f(-z-z_0)$, and $-\overline{f(\bar z-z_0)}=-\bar f(z-\bar z_0)$. The holomorphic field (call it $F(z)$) corresponding to the full potential $\phi$ is then just the sum of these individual 'reflected' holomorphic copies of $f(z-z_0)$. To satisfy the boundary conditions specified, we need to ensure that the sum of all $f$ and its reflections is purely imaginary plus a constant real part on any of the circular boundaries of electrodes. Because of the extra symmetry brought using the conformal mapping, we can focus on just one of these boundary circles. 
For simplicity, let's choose the factor $A$ so that the boundary circles all have radius 1. Assuming the (locally) holomorphic function $f(z)$ has a singular expansion about $z=0$ (i.e. $f(z-z_0)$ is [potentially] singular near $z=z_0$), we can write $f(z)=\alpha \ln z + \sum_{n\in \mathbb Z} c_n z^n$. I'll also assume that $c_n=0$ for $n> 0$.
Calling $g(z-z_0)$ the sum of the three reflections of $f(z-z_0)$, we assume $g(z)$ is analytic at $z=0$, so we can write $g(z)=\sum_{n\geq 0}b_n z^n$. In order for the real part of $f(z)+g(z)$ to be constant on the unit circle, we need each nonzero frequency sector in the "Fourier expansion" to be purely imaginary. Since $f(z)$ has a singular expansion while $g(z)$ has an analytic expansion, this means imposing conditions on $b_n$ and $c_{-n}$. In order for $b_nz^n+c_{-n}z^{-n}\big|_{\partial \mathbb D}=b_ne^{in\theta}+c_{-n}e^{-in\theta}$ to be imaginary, we need $c_{-n}=-\bar b_n=-\frac{1}{n!}\overline{g^{(n)}(0)}$, where $g^{(n)}$ is the $n$-th complex derivative of $g$. 
From this, we can write 
\begin{align*}
f(z)=\alpha\ln z-\sum_{n\geq 0}\frac{1}{n!}\overline{g^{(n)}(0)}z^{-n},
\end{align*}
and since $\overline{g(z)}=\sum_{n=0}^\infty \frac{1}{n!}\overline{g^{(n)}(0)}\overline{z^n}$, we have $f(z)=\alpha\ln z-\overline{g(\bar z^{-1})}=\alpha\ln z - \bar g(z^{-1})$. Using the original definition of $g(z-z_0)=-\bar f(z-\bar z_0)-f(-z-z_0)+\bar f(-z-\bar z_0)$, we have a closed (functional) equation for $f(z)$:
\begin{align*}
f(z)=\alpha\ln z + f(\bar z_0-z_0+z^{-1})+\bar f(-2\bar z_0-z^{-1})-f(-\bar z_0-z_0-z^{-1})
\end{align*}
In cases where the radii of the electrodes is very small, it is reasonable to keep only the first term $\alpha \ln z$, so that $\psi(x,y)=\alpha\ln|z|$. Then  $\phi(x,y)$ is given by $\alpha (\ln |z-z_0|+\ln|z+\bar z_0|-\ln|z+z_0|-\ln|z-\bar z_0|)$, and the potential in the original coordinates ($s$ say) is given by $\phi(\sigma(s))$, where $\sigma$ was the Moebius transformation defined earlier. This approximation for $\phi$ (in the original domain) might also be a good starting point for numerical relaxation methods.
For finite radius electrodes, it is possible to use the above relation for $f$ iteratively to improve the accuracy of the approximation. The expansion can be represented diagrammatically as a sum over nodes of a sort of truncated tree, where each leaf is associated with a logarithm of a finite continued fraction. 
