n = 38;

nsq = n*n;

u=1:n;
v=u;

for k=1:n
  topboundary(k)=exp(-u(k)*v(1));
  leftboundary(k)=exp(-u(1)*v(k));
  rightboundary(k)=exp(-u(k)*v(n));
  bottomboundary(k)=exp(-u(n)*v(k));
end

a = zeros(nsq,nsq);

for k=1:nsq
  a(k,k)=4;
end
for k=1:nsq-1
  if ((k/n)~=floor(k/n)) a(k,k+1)=-1; end;
end
for k=2:nsq
  if (((k-1)/n)~=floor((k-1)/n)) a(k,k-1)=-1; end;
end

for k=1:nsq-n
    a(k,k+n) = -1;
end

for k=n+1:nsq
    a(k,k-n) = -1;
end

b = zeros(nsq,1);

for k=2:n
  b(k) = topboundary(k);
  b(n*(k-1)+1) = leftboundary(k);
  b(n*(k-1)) = rightboundary(k);
  b(nsq-n+k) = bottomboundary(k);
end

b(1) = topboundary(1)+leftboundary(1);
b(n) = topboundary(n)+rightboundary(1);
b(nsq-n+1) = bottomboundary(1)+leftboundary(n);
b(nsq) = bottomboundary(n)+rightboundary(n);

x = a\b;

for i=1:n
  for j=1:n
     heat(i,j)=x(i+(j-1)*n);
  end
end

contour(u,v,heat,40);
