function u_val = heat(pos,rectangle,T,u0,f,bdry,h);

% u_val = heat(pos,rectangle,T,u0,f,bdry,h)
%
% solves the heat equation u_t - \Delta u(x) = f with dirichlet 
% boundary conditions on rectangle [-a,a] x [-b,b] using a five-point 
% stencil and explicit Euler time-stepping
%
%   pos         [x,y] a point of the rectangle
%   rectangle   [a,b]
%               2a   length in x-direction
%               2b   length in y-direction
%   T           final time
%   u0          intial value                        (constant)
%   f           right-hand side                     (constant)
%   bdry        Dirichlet data [xl,xr,yl,yu]
%               xl  boundary value on {-a} x [-b,b]   (constant)
%               xr  boundary value on { a} x [-b,b]   (constant)
%               yl  boundary value on [-a,a] x {-b}   (constant)
%               yu  boundary value on [-a,a] x { b}   (constant)
%   h   h = 1/n suggested approximate grid-size
%   u_val       solution at pos = [x,y]

% the grid
a = rectangle(1); b = rectangle(2);
n = 2*ceil(1/2/h); % make n even
nx = ceil(2*a)*n+1; ny = ceil(2*b)*n+1;
dx = 2*a/(nx-1); dy = 2*b/(ny-1);
nt = ceil(4*T)*ceil(1/h^2); dt = T/nt; 
x = 1:nx; y=1:ny;

% initial and boundary values
u = u0*ones(nx,ny); 
u(x(1),y) = bdry(1); u(x(end),y) = bdry(2);
u(x,y(1)) = bdry(3); u(x,y(end)) = bdry(4);

% the five-point stencil
stencil = [-1/dy^2 -1/dx^2 2*(1/dx^2+1/dy^2) -1/dx^2 -1/dy^2];

% the time-stepping
x = x(2:end-1); y = y(2:end-1); 
for k=1:nt
    u(x,y) = u(x,y) + dt*(f - stencil(1)*u(x,y-1)...
                            - stencil(2)*u(x-1,y)...
                            - stencil(3)*u(x,y)...
                            - stencil(4)*u(x+1,y)...
                            - stencil(5)*u(x,y+1));
end

% the solution
u_val = u(1+round((pos(1)+a)/dx),1+round((pos(2)+b)/dy));

return;