function vnew=dampj(vold,f)
% DAMPJ Damped Jacobi code for Helmholtz
%
% Tim's damped Jacobi code
% Check out the fancy matlab!
%
global sigma
%
% Stay awake! The range of unknowns is 2 to nl-1. MATLAB does not
% support zero-based arrays.
%
omega=2/3;
nl=length(vold); 
h=1/(nl-1); h2=h*h;
d=omega/(2 + h*h*sigma);
vnew=vold;
%
% the matlab vector instructions do what the loop does, only faster
%
%for j=2:nl-1
%	vnew(j)=vold(j-1)+vold(j+1)+h2*f(j);
%end
%
vnew(2:nl-1)=(1-omega)*vold(2:nl-1) ...
             + d*(vold(1:nl-2)+vold(3:nl)+h2*f(2:nl-1));