function Heat2DStatic()
nx = 200; ny = 200;
%nx = 2*500; ny = nx;
function res = f(x)
  res = x(:,1).^2;
end%function

%%%%%%%%%%%%%%%%%%%%%%%%%%
hx = 1/(nx+1); hy = 1/(ny+1);
Dxx = spdiags([ones(nx,1) -2*ones(nx,1) ones(nx,1)],[-1 0 1],nx,nx)/(hx^2);
Dyy = spdiags([ones(ny,1) -2*ones(ny,1) ones(ny,1)],[-1 0 1],ny,ny)/(hy^2);
A = -kron(Dxx,speye(ny))-kron(speye(nx),Dyy);
x = hx:hx:1-hx; y = hy:hy:1-hy; 
[xx,yy] = meshgrid(x,y);
fvec = f([xx(:),yy(:)]);
t0 = cputime();
u = A\fvec;
solutionTime = cputime()-t0

%%% add on the zero boundary values for a nicer graphics
x = 0:hx:1; y = 0:hy:1;
[xx,yy] = meshgrid(x,y);  uu = zeros(size(xx));
uu(2:ny+1,2:nx+1) = reshape(u,ny,nx);
mesh(xx,yy,uu)
xlabel('x'); ylabel('y');
end%function
