n=20; % numer of r-divisions
m=80; % number of theta-divisions
% Note: going all the way to the boundary will result in
% divisions by zero.

rho=1; around=2*pi;
dr=rho/(n+1); dth=around/m; dfi=around/m;
r=0:dr:rho; th=0:dth:around;

[rg,thg]=meshgrid(r,th);

% This is a very simple quadrature rule to do the integration.
% Simpson or something else would be just as easy to implement.
% I leave it as a exercise. :o)
usol=0; %% I'm not sure if the for loop is necessary...
for k=0:m-1,
  usol=usol+(rho^2-rg.^2).*fu(k*dfi)./(rho^2-2*rho*rg.*cos(thg-k*dfi)+rg.^2);
end
usol=usol/m;

X=rg.*cos(thg); 
Y=rg.*sin(thg);

surf(X,Y,usol);
rotate3d on;
colormap summer; shading faceted;
%axis ij;
xlabel('x');ylabel('y');
