function sol=helmholtz(f,repeats);
% HELMHOLTZ Solve discrete Helmholtz equation with zero boundary data
%
global sigma 
nt=length(f); nl=length(f)-1; nv=nl-1;
rhs=f(2:nt-1);
sol=zeros(nt,1);
%
% a is the coefficent matrix of the differential equation
% this lets us compute the exact solution for any right side
% and make sure that everything else is ok
%
% I build a with the matlab sparse matrix commands.
%
hm2=nl*nl;
da=(sigma+2*hm2)*ones(nv,1); od=hm2*ones(nv,1);
a=spdiags([-od, da, -od], [-1,0,1], nv,nv);
air=a\rhs;
sol(2:nt-1)=air;
if nargin==2
   for i=2:repeats
       air=a\rhs;
   end
end