%%%%% script file to solve the heat equation on a unit square
clear *
n = 4*100;
f = @(x,y) x.^2; % describe the heating contribution
%%%%%%%%%%%%%%%%%%%%%%%%%%
%%no modifications necessary beyond this line
h = 1/(n+1);
Dxx = spdiags(ones(n,1)*[-1 2 -1],[-1 0 1],n,n)/h^2;
A = kron(Dxx,eye(n)) + kron(eye(n),Dxx);
x = h:h:1-h; y = x;
[xx,yy] = meshgrid(x,y);
fvec = f(xx(:),yy(:));
tic(); t0 = cputime();
u = A\fvec;
for jj = 1:20
  fVec = A*u;
  u = A\fvec;
endfor
toc()
cputime - t0
