function Eg7_3Rev()
% Simulates the cooling of a plate whose boundary temperature is fixed
% at zero.

% Establish the Initial Temperature Distribution...
a = 0; b = 6; n = 61; x = linspace(a,b,n);
c = 0; d = 4; m = 41; y = linspace(c,d,m);
TVals = fOnGrid(x,y,@T);
TVals(1,:) = zeros(1,n); TVals(m,:) = zeros(1,n);
TVals(:,1) = zeros(m,1); TVals(:,n) = zeros(m,1);
nSteps = 200;
close all
figure
colormap('hot')
for tau=0:3*nSteps
% Display the current temperature distribution...
if rem(tau,10)==0
pcolor(TVals)
shading interp
caxis('manual')
title(sprintf('tau = %1d',tau),'Fontsize',14)
axis equal off
pause(.1)
shg
end
% Update the temperature distribution...
TVals = Average(TVals);
end

function tau = T(x,y)
tau = 100*exp(-.4*( (x-1)^2 + 0.7*(y-3)^2)) + ...
80*exp(-.2*(2*(x-5)^2 + 1.5*(y-1)^2));

function fVals = fOnGrid(x,y,f)
% x is a 1-by-n vector
% y is a 1-by-m vector
% f is a function handle that identifies a function of two variables.
% fVals is an m-by-n matrix where fVals(i,j) = f(x(j),y(i)).

n = length(x); m = length(y); fVals = zeros(m,n);
for j=1:n
for i = 1:m
fVals(i,j) = f(x(j),y(i));
end
end

function B = Average(A)
% A is an m-by-n matrix
% B is an m-by-n matrix
% The first and last row and first and last column of B is zero.
% Otherwise, B(i,j) is the average of A(i,j)'s four neighbors, i.e.,
% A(i+1,j), A(i,j+1), A(i-1,j), and A(i,j-1).

[m,n] = size(A);
B = zeros(m,n);
for i=2:m-1
for j=2:n-1
North = A(i-1,j);
East = A(i,j+1);
South = A(i+1,j);
West = A(i,j-1);
B(i,j) = (North+East+South+West)/4;
end
end