clear;

r1=0; dr=0.05; r2=1;
a1=0; da=2*pi/100; a2=2*pi;

rint=r1:dr:r2;
aint=a1:da:a2;

% Three pointlike heat sources on the unit disc boundary
dev=0.1;
for i=1:3
  a(i)=sin(i*2*pi/3+dev);b(i)=cos(i*2*pi/3+dev);
end

X=rint'*cos(aint); 
Y=rint'*sin(aint); 

pk=0;
for i=1:3
  num=1/(2*pi)*(1-X.^2-Y.^2);
  denom=((X-a(i)).^2+(Y-b(i)).^2);
  pk=pk+num./denom;
end

surf(X,Y,pk); shading interp;
axis([-1 1 -1 1 0 5])
colormap hot
rotate3d
