
%well location x,y coordinates
sx = [-1 0 1]; sy = [1 2 1];

%observation locations
px = [-1 0 1];  py = [0 0 0];

%solve inverse problem

U = [.0680 .1106 .1214]';  % observation data

for i = 1:3, for j=1:3,
    G(i,j) = (1/(2*pi))*( besselk(0, sqrt( (px(i)-sx(j))^2 + (py(i)-sy(j))^2)) + ...
               besselk(0, sqrt( (px(i)-sx(j))^2 + (py(i)+sy(j))^2)) );
end; end;

sigma = G\U,

