hr = 1/Nr; ht = 1/Nt;
tau = input('tau= ');
lambda_r = tau/hr^2; lambda_t = tau/ht^2;
I = eye(Nt);
J = diag(ones(Nt-1,1),1)+diag(ones(Nt-1,1),-1)+diag(1,Nt-1)+diag(1,-Nt+1);
r = 1/Nr:1/Nr:(Nr-1)/Nr;
di = 1 + 2*lambda_t./(r.*r);
ti = -lambda_t./(r.*r);
a = diag(di); b = diag(ti);
B = kron(I,a) + kron(J,b);
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;

fi = 1-2*lambda_r;
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;
d = fi * eye(Nr-1) + diag(ri,1) + diag(li,-1);
C = kron(I,d);
D = zeros((Nr-1)*Nt+1); D(2:(Nr-1)*Nt+1,2:(Nr-1)*Nt+1) = C;
g = zeros(1,(Nr-1)*Nt+1); g(1,1) = 1-4*lambda_r;
for j = 0:Nt-1;
    g(1,(Nr-1)*j+2) = 4*lambda_r/Nt;
end;
D(1,1:(Nr-1)*Nt+1) = g;
F = inv(A)*D;
opts.disp = 0;
abs(eigs(F,1,'lm',opts))