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;