clear;
N=100;
MB=8*N/10;
MT=9*N/10;
ML=4*N/10;
MR=9*N/10;

CB=2*N/10;
CT=4*N/10;
CL=2*N/10;
CR=3*N/10;
maxChange=1000;
T(1:N,1:N)=0;
newT(1:N,1:N)=0;

T(MB:MT,ML:MR)=400;
T(CB:CT,CL:CR)=1200;
T(1,1:N)=80;
T(N,1:N)=80;
T(1:N,1)=80;
T(1:N,N)=80;


while maxChange>1

	for i=2:N-1
	  for j=2:N-1
	    newT(i,j)=(T(i-1,j)+T(i+1,j)+T(i,j+1)+T(i,j-1))/4;
	  end
	end
	newT(MB:MT,ML:MR)=400;
	newT(CB:CT,CL:CR)=1200;
	newT(1,1:N)=80;
	newT(N,1:N)=80;
	newT(1:N,1)=80;
	newT(1:N,N)=80;
	maxChange=sum(sum(abs(T-newT)));
	T=newT;

end
mesh(1:N, 1:N, T)

