U=zeros(n+1,n+1);
Unew=zeros(n+1,n+1);
U(2:n,2:n)=randn(n-1,n-1);
surf(U),pause
for iterations = 1:50,
  for j=2:n,
    for k=2:n,
Unew(j,k) = (1-theta)*U(j,k)+(theta/4)*(U(j+1,k)+U(j-1,k)+U(j,k+1)+U(j,k-1));
    end;
  end;
  surf(Unew);U=Unew;pause,
end