hr=1/Nr;ht = 1/Nt;
tau = input(' tau= ');
lambda_r=tau/hr^2;lambda_t=tau/ht^2;
r=hr:hr:(Nr-1)/Nr;
s=1/Nr:1/Nr:(Nr-2)/Nr;
ri = lambda_r + (lambda_r*hr)./2*s;
u=2/Nr:1/Nr:(Nr-1)/Nr;
li = lambda_r -(lambda_r*hr)./2*u;
di=1-2*lambda_r-(2*lambda_t./(r.*r));
ti = lambda_t./(r.*r);
a = diag(di)+diag(ri,1) + diag(li,-1);
b = diag(ti);
J =diag(ones((Nt-1),1),1) + diag(ones((Nt-1),1),-1) 
            +diag(1,Nt-1)+diag(1,-(Nt-1));
I = eye(Nt);
B = kron(J,b) + kron(I,a);
A = zeros((Nr-1)*Nt+1); A(2:(Nr-1)*Nt+1,2:(Nr-1)*Nt+1) = B;
c = zeros(1,(Nr-1)*Nt+1);
c(1,1) = 1 - 4 * lambda_r;
for i = 0:Nt-1;
      c(1,(Nr-1)*i+2) = 4 * lambda_r/Nt;
end;
A(1,1:(Nr-1)*Nt+1) = c;
opts.disp = 0;
abs(eigs(A,1,'lm',opts))