Ilustrace metody fiktivních nábojů pro dva vodiče ve tvaru koule

Příklad se počítá na cvičení z elektrodynamiky 
(viz http://utf.mff.cuni.cz/~ledvinka/?278656 a http://utf.mff.cuni.cz/~ledvinka/PrElektrodynamika/Kapacity.pdf)

(* funkce počítají to, co my na cvičení, jen kvůli pohodlí začínám s indexy od jedné *)
spoctiQS[a_,b_,d_,U1_,U2_,n_Integer]:=Block[{s,S,q,Q},
s[1]=0;
S[1]=0;
q[1]=a U1;
Q[1]=b U2;
Do[
s[i]=a^2/(d-S[i-1]);
S[i]=b^2/(d-s[i-1]);
q[i]= -a/(d-S[i-1]) Q[i-1];
Q[i]= -b/(d-s[i-1]) q[i-1];
,{i,2,n}];
Transpose[ Table[{q[i],Q[i],s[i],S[i]},{i,1,n}]]
]
Phi[a_,b_,d_,U1_,U2_,n_Integer]:=Block[{s,S,q,Q,i,tmp},
{q,Q,s,S}=spoctiQS[a,b,d,U1,U2,n];
tmp=Evaluate[ Sum[ q[[i]]/Sqrt[x^2+(z-s[[i]])^2]+Q[[i]]/Sqrt[x^2+(z-d+S[[i]])^2],{i,1,n}] ];
Function[{x,z},Evaluate[tmp]]
]
maticeKapacit[a_,b_,d_,n_Integer]:=Block[{s,S,q,Q,U1,U2,i,tmp},
{q,Q,s,S}=spoctiQS[a,b,d,U1,U2,n];
D[{Total[q],Total[Q]},{{U1,U2}}]
]

maticeKapacit[a,b,d,4]//Simplify//MatrixForm (*Bez Simplify není vidět, že (Subscript[C, ab])je symetrická*)

vzdalenost=3;
polomer1=1;
polomer2=1;

ϕ=Phi[polomer1,polomer2,vzdalenost,-1,1,5];

ContourPlot[ ϕ[x,z],{x,-2,2},{z,-2,2+vzdalenost},AspectRatio->Automatic,PlotPoints->20,
Contours->40,PlotRange->{-1,1},
RegionFunction->Function[{x,z},x^2+z^2>polomer1^2 &&x^2+(z-vzdalenost)^2>polomer2^2 ],BoundaryStyle->{Thick, Black}]

(*Ukázka vlivu počtu fiktivních nábojů*)
Table[ ContourPlot[ Phi[polomer1,polomer2,vzdalenost,-1,1,n][x,z],{x,-2,2},{z,-2,2+vzdalenost},AspectRatio->Automatic,PlotPoints->20,
Contours->20,ImageSize->300,PlotRange->{-1,1},PlotLabel->"n="<>ToString[n],
RegionFunction->Function[{x,z},x^2+z^2>polomer1^2 &&x^2+(z-vzdalenost)^2>polomer2^2 ],BoundaryStyle->{Thick, Black}],
{n,1,4}]//GraphicsRow

Varianta s U1=0, U2=1

Table[ ContourPlot[ Phi[polomer1,polomer2,vzdalenost,0,1,n][x,z],{x,-2,2},{z,-2,2+vzdalenost},AspectRatio->Automatic,PlotPoints->20,
Contours->20,ImageSize->300,PlotRange->{0,1},PlotLabel->"n="<>ToString[n],
RegionFunction->Function[{x,z},x^2+z^2>polomer1^2 &&x^2+(z-vzdalenost)^2>polomer2^2 ],BoundaryStyle->{Thick, Black}],
{n,1,4}]//GraphicsRow