with(plots):
with(plottools):
com:=(x,y)->x+I*y;
m:=2;
Z[1]:=com(1,1);
Z[2]:=com(3,2);
q[1]:=3;
q[2]:=5;
phi := (z) -> sum(q[i]*log(z-Z[i]), i=1..m);
# contourplot(Re(phi(com(x,y))),x=0..4,y=0..4,grid=[50,50]);
# contourplot(Re(phi(com(x,y))),x=-10..15,y=-10..15,grid=[70,70],thickness=2,coloring=[red,blue]);
# contourplot(Re(phi(com(x,y))),x=-10..15,y=-10..15,thickness=2,coloring=[red,blue]);

me_init := proc(pmax,a,Q)
    local k;
    for k to pmax do
	a[k] := -sum(q[i]*Z[i]^k/k, i=1..m);
    od;
    Q := sum(q[i], i=1..m);
end;

me := (p,z) -> Q*log(z) + sum(a[k]/z^k, k=1..p);

me_init(10,a,Q);

gr:=160;
contourplot(Re(phi(com(x,y))),x=-15..15,y=-15..15,thickness=2,coloring=[blue,red],contours=[seq(3*i,i=-10..10)],grid=[gr,gr],scaling=CONSTRAINED,title=`exact potential`);
contourplot(Re(me(0,com(x,y))),x=-15..15,y=-15..15,thickness=2,coloring=[blue,red],contours=[seq(3*i,i=-10..10)],grid=[gr,gr],scaling=CONSTRAINED,title=`multipole, p=0`);
contourplot(Re(me(1,com(x,y))),x=-15..15,y=-15..15,thickness=2,coloring=[blue,red],contours=[seq(3*i,i=-10..10)],grid=[gr,gr],scaling=CONSTRAINED,title=`multipole, p=1`);
contourplot(Re(me(2,com(x,y))),x=-15..15,y=-15..15,thickness=2,coloring=[blue,red],contours=[seq(3*i,i=-10..10)],grid=[gr,gr],scaling=CONSTRAINED,title=`multipole, p=2`);
contourplot(Re(me(3,com(x,y))),x=-15..15,y=-15..15,thickness=2,coloring=[blue,red],contours=[seq(3*i,i=-10..10)],grid=[gr,gr],scaling=CONSTRAINED,title=`multipole, p=3`);
contourplot(Re(me(10,com(x,y))),x=-15..15,y=-15..15,thickness=2,coloring=[blue,red],contours=[seq(3*i,i=-10..10)],grid=[gr,gr],scaling=CONSTRAINED,title=`multipole, p=10`);