function PlateDynamic()
%%%%% function file to solve the dynamic heat equation on a unit square
T = 0.5;  dt = 0.02;

nx = 34; ny = 35;
function res = f(x)
  res = x(:,1).^2;
end%function

%%%%%%%%%%%%%%%%%%%%%%%%%%
t = 0:dt:T;   utrace = zeros(size(t));

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);
A = speye(nx*ny,nx*ny) + dt*A;
[R,p,P] = chol(A); % Cholesky factorization
% R = chol(A);
x = hx:hx:1-hx; y = hy:hy:1-hy; 
u = zeros(nx*ny,1);  % define the initial temperature
[xx,yy] = meshgrid(x,y);  fvec = dt*f([xx(:),yy(:)]);

tic
for k = 2:length(t);
  u = P*(R\(R'\(P'*(u+fvec)))); % one time step
%  u = R\(R'\(u+fvec)); % one time step
  utrace(k) = u(55);
end
solutionTime = toc

%%% 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);

set (0,'defaultlinelinewidth',2) 

figure(1); surf(xx,yy,uu); xlabel('x'); ylabel('y');
figure(2); plot(t,utrace); grid on ;xlabel('Time t'); ylabel('Temp u');

end%function
