L'Hospital came up with separation of variables in 1750, and it is now the physicist's handiest tool for solving partial differential equations. Using Laplace's equation and given the appropriate boundary conditions, it is possible to determine V (the electric potential) or W (the scalar potential of the auxiliary field) using separation of variables. In general, given a partial differential equation in a function Φ(x_{1}, x_{2}, ..., x_{n}), we break it up into a set of independent ordinary differential equations such that

Φ(x_{1}, x_{2}, ..., x_{n})=X_{1}(x_{1})X_{2}(x_{2})...X_{n}(x_{n})

X_{i} can be easily solved and substituted back in. This is not possible for many Φ, but for lots of problems in physics this process is capable of finding solutions from specified boundary conditions. The solutions for V are unique for a given set of boundary conditions, so the answer we get through this process is the actual solution.

It is important to pick the right coordinate system. Let's pick Cartesian coordinates.

V(x,y,z)=X(x)Y(y)Z(z)

∂^{2}V/∂x^{2}+∂^{2}V/∂y^{2}+∂^{2}V/∂z^{2}=0 *Laplace's Equation*

∂^{2}X/∂x^{2}*Y(y)Z(z)+∂^{2}Y/∂y^{2}*X(x)Z(z)+∂^{2}Z/∂z^{2}*X(x)Y(y)=0 *Plugging V into Laplace's Equation*

∂^{2}X/∂x^{2}*1/X(x)
+∂^{2}Y/∂y^{2}*1/Y(y)+∂^{2}Z/∂z^{2}*1/Z(z)=0 *Dividing by V*

∂^{2}X/∂x^{2}*1/X=C_{1}, ∂^{2}Y/∂y^{2}*1/Y=C_{2}, ∂^{2}Z/∂z^{2}*1/Z=C_{3} C_{1}+C_{2}+C_{3}=0 *Follows from above*

Depending on your charge distribution, you can pick C_{1}, C_{2} and C_{3} however you like. Let's assume that the potential goes to 0 as x goes to infinity. Then we will pick C_{1} to be positive and C_{2}, C_{3} to be negative.

C_{1}=k^{2}+l^{2}, C_{2}=-k^{2}, C_{3}=-l^{2}

∂^{2}X/∂x^{2}=(k^{2}+l^{2})X, ∂^{2}Y/∂y^{2}=-k^{2}Y,
∂^{2}Z/∂z^{2}=-l^{2}Z,

Now we have our ordinary differential equations. Solving...

X(x)=Ae^{√(k2+l2) x}+Be^{-√(k2+l2) x}

Y(y)=C sin(ky)+ D cos(ky)

Z(z)=E sin(lz)+ F cos(lz)

While the formula V(x,y,z)=(Ae^{√(k2+l2) x}+Be^{-√(k2+l2) x})(C sin(ky)+ D cos(ky))(E sin(lz)+ F cos(lz)) might not seem particularly enlightening, with the appropriate boundary conditions you can solve for A, B, C, D, E, F, k and l, many of which turn out to be zero. It might also seem like the sine and cosine terms give infinite solutions, making this whole process worthless; but it is exactly these infinite solutions that allows us to find a solution for V.

Laplace's Equation is linear, such that any linear combination of solutions to Laplace's equation is also a solution to Laplace's Equation. Thus our solution ends up being an infinite sum ∑_{n, 0->∞} V_{n}, where k and l are linear functions of n. The completeness and orthogonality of the sinusoidal functions allows us to use this infinite sum to find the actual V, not unlike a Fourier series. In practice, most of the terms usually drop out or the contributions of successive terms decreases quickly such that you can get a reasonable approximation for V with only a few terms.

This process also works in a whole bunch of other coordinate systems. The general solution in spherical coordinates (assuming azimuthal symmetry) is:

V(r, θ)=∑_{n, 0->∞}(A_{n}r^{n}+B_{n}/r^{n+1})P_{n}(cos θ)

Where P_{n}(x) are the Legendre polynomials.

Sources:

http://mathworld.wolfram.com/SeparationofVariables.html

Griffiths, D. *Introduction to Electrodynamics* 3rd edition © 1999