function H = jacobi(s)
% JACOBI : Demonstrate Jacobi iteration on 2D temperature problem
% 
% With no arguments, jacobi pauses after each iteration
% H = jacobi(s) does s iterations

m = 100;   % size of interior of grid
if nargin == 0
    pausing = true;
    s = 1000000;
else
    pausing = false;
end;

% Initialize 2d temperature array including boundary conditions
H = zeros(m+2,m+2);               % 0 degrees C (32 F) at walls
radiator = round(m/2 - .2*m) : round(m/2 + .2*m);
H(1,radiator) = 100;              % 100 degrees C (212 F) at radiator
if pausing
    surfc(H), shg;
    pause;
end;

% Indices for interior points and neighbor points
this = 2:m+1;
prev = this-1;
next = this+1;

% Now the main loop
for step = 1:s
    
    H(this,this) = (H(this,prev) + H(prev,this) + H(this,next) + H(next,this))/4;
    if pausing
        step
        surfc(H),shg;
        pause;
    end;
end;

surfc(H),shg