function v = bem_solve(xx,h,x0,k,N)
%
% function v = bem_solve(xx,h,x0,k,N)
%
% This function solves the boundary integral equation (2.33) in the notes by the simplest
% boundary element method in the 2D case (though the code applies in 3D if
% the function G is redefined to be the 3D Green's function).
%
% The method used is that explained in the last problem on the problem
% sheet.
%
% On entry:
%
% xx    is an N x 2 matrix, the jth row of xx containing the position
%       vector of the collocation point in the jth element
% h     is an N x 1 vector, the jth entry of h the length of the jth
%       element
% x0    is a 1x2 vector, the position vector of the point source of sound
% k     is the wavenumber
% N     is the number of boundary elements
%
% On exit:
%
% v     is an N x 1 (complex) vector, v(j) an approximation to the normal
%       derivative of the pressure (in the direction pointing into the domain) on
%       element j
%

% First compute the matrix and right hand side of the linear system to be
% solved.
%
A = zeros(N,N);
b = zeros(N,1);
for j = 1:N
    b(j) = G(xx(j,:),x0,k);
end
for i = 1:N
    for j = 1:N
        if i~=j
            A(i,j) = h(j)*G(xx(j,:),xx(i,:),k);
        end
    end
end
%
% Now solve the linear system
%
v = -A\b;