Laplace's Equation in 3 dimensions is given by $$\nabla^2f=\frac{ \partial^2f}{\partial x^2}+\frac{ \partial^2f}{\partial z^2}+\frac{ \partial^2f}{\partial y^2}=0$$ and is a very important PDE in many areas of science. One of the usual ways to solve it is by seperation of variables. We let $f(x,y,z)=X(x)Y(y)Z(z)$ and then the PDE reduces to 3 independent ODE's of the form $X''(x)+k^2 x=0$ with $k$ a constant.
This method works for a surprising number of coordinate systems. I can use cylindrical coordinates, spherical coordinates, bi-spherical coordinates and more. In fact in 2 dimensions, I can take $\mathbb{R}^2$ to be $\mathbb{C}$, and then any analytic function $f(z)$ maps the Cartesian coordinates of $\mathbb{R}^2$ onto a set of coordinates that is separable in Laplace's equation. This follows from analytic functions also being harmonic functions. For example $f(z)=z^2$ maps the cartesian coordinates to the parabolic coordinates.
So how about in 3 dimensions? Is there a way to describe all coordinate systems such that the Laplacian separates in this way?I have no idea how to start. I think that Confocal Ellipsoidal Coordinates will work but I'm not sure how to verify it. I'm not familIar with differential geometry, so any help appreciated.