Dave, your method is a version of Smythe's "oblate spheroidal harmonics"
method where he, too, makes a transformation to boundary matching
coordinates (his section 5.271)
His handling of the disc (pg 111-114, section 5.00 & ff) is quite
different and would interest you. Because you'll probably have difficulty
finding Smythe, let me try to outline at least the beginning of his
method:
The equation of a surface is F(x,y,z) = C; for each value of C this will
be used as an equipotential surface: V = f(C).
He then applies Laplace's equation DEL^2 V = 0 and gets
DEL^2 V = f''(C)*{Grad C}^2 + f'(C)*DEL^2 C=0 =>
DEL^2 C/{Grad C}^2 = - f''(C)/f'(C) = G(C) That this is a function only
of C is then the requirement that F(x,y,z)=C can be an equipotential .
He then integrates this to get the potential (much skipped):
V = f(C) = A INT exp{-INT G(C) dC} dC + B. A & B are determined by
specifying the potential on two of the surfaces.
He then applies this to the "nonintersecting confocal conicoids"
x^2/(a^2+C) + y^2/(b^2+C) + z^2/c^2+C) = 1 and explicitly solves the disc
as a special case.
In his section 5.271 he develops the "oblate spheroidal coordinates" by
taking
b=c , y = r*cos(phi) and z = r* sin(phi) . . . .
I hope you can get Smythe, cuz I'm sure the above is over-simplified and
bungled. (I could mail you a photocopy, if needed.)