syms r1 r2 x y z a b c1 c2 Vzx Vzy Vzz Vxx Vxy Vyy Bx By Bz Mz Mx My u0 f1 f2 f3 f4 f5 f6 f7 f8 Bxx Bxy Bxz Byy Byz;
r1=sqrt((x-a)^2+(y-b)^2+(z-c1)^2);
r2=sqrt((x-a)^2+(y-b)^2+(z-c2)^2);
Vzx=(x-a)*(1/r2^3-1/r1^3);
Vzy=(y-b)*(1/r2^3-1/r1^3);
Vzz=(c2-z)/r2^3-(c1-z)/r1^3;
Vxx=((y-b)^2-(x-a)^2)/((y-b)^2+(x-a)^2)^2*((c1-z)/r1-(c2-z)/r2)-(x-a)^2/((y-b)^2+(x-a)^2)*((c1-z)/r2^3-(c2-z)/r1^3);
Vxy=2*(x-a)*(y-b)/((y-b)^2+(x-a)^2)^2*((z-c1)/r1-(z-c2)/r2)-(x-a)*(y-b)/((y-b)^2+(x-a)^2)*((z-c1)/r2^3-(z-c2)/r1^3);
Vyy=((x-a)^2-(y-b)^2)/((y-b)^2+(x-a)^2)^2*((c1-z)/r1-(c2-z)/r2)-(y-b)^2/((y-b)^2+(x-a)^2)*((c1-z)/r2^3-(c2-z)/r1^3);
f1=Bx-u0/(4*pi)*(Mx*Vxx+My*Vxy+Mz*Vzx);
f2=By-u0/(4*pi)*(Mx*Vxy+My*Vyy+Mz*Vzy);
f3=Bz-u0/(4*pi)*(Mx*Vzx+My*Vzy+Mz*Vzz);
f4=Bxx-diff(Bx,x);
f5=Bxy-diff(Bx,y);
f6=Bxz--diff(Bx,z);
f7=Byy-diff(By,y);
S=solve(f1==0,f2==0,f3==0,f4==0,f5==0,f6==0,f7==0,a,b,c1,c2,Mx,My,Mz);