% Variables related to the numerical integration procedure (Simpson's rule)
% The interval [0, 2 pi] is divided into 200 subintervals:
hvec = 2*pi/600 * [1 repmat([4 2], 1, 99) 4 1]';
omega = 0:2*pi/200:2*pi;
cx = cos(omega);
cy = sin(omega);

% Discretization of the domain of the solution
nr = 20;
nt = 100;
u = zeros(nr, nt+1);
r = 0: 1/nr :(1-1/nr);
theta = 0: 2*pi/nt :2*pi;
X = r' * cos(theta);
Y = r' * sin(theta);

for ir = 1:nr
  for it = 1:(nt+1)
    num = (1 - X(ir, it)^2 - Y(ir, it)^2) .* f(omega);
    denom = 2*pi*((X(ir, it) - cx).^2 + (Y(ir, it) - cy).^2);
    u(ir, it) = (num ./ denom) * hvec;
  end
end

surf(X,Y,u);

