hr = 1/Nr; ht = 1/Nt;
tau = input('tau= ');
lambda_r = tau/hr^2; lambda_t = tau/ht^2;
lambda1_r=lambda_r/2; lambda1_t = lambda_t/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*lambda1_t./(r.*r);
ti = -lambda1_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*lambda1_r;
for i = 0:Nt-1;
    c(1,(Nr-1)*i+2) = 4*lambda1_r/Nt;
end;
A(1,1:(Nr-1)*Nt+1) = c;

fi = 1-2*lambda1_r;
s = 1/Nr:1/Nr:(Nr-2)/Nr;
ri = lambda1_r + (lambda1_r*hr)./2*s;
u = 2/Nr:1/Nr:(Nr-1)/Nr;
li = lambda1_r - (lambda1_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*lambda1_r;
for j = 0:Nt-1;
    g(1,(Nr-1)*j+2) = 4*lambda1_r/Nt;
end;
D(1,1:(Nr-1)*Nt+1) = g;
F = inv(A)*D;

di = 1 + 2*lambda1_r;
s = 1/Nr:1/Nr:(Nr-2)/Nr;
ri = -lambda1_r - lambda1_r * hr./2*s;
u = 2/Nr:1/Nr:(Nr-1)/Nr;
li = -lambda1_r + lambda1_r * hr./2*u;
a = di*eye(Nr-1) + diag(ri,1) + diag(li,-1);
G = kron(I,a);
H = zeros((Nr-1)*Nt+1); H(2:(Nr-1)*Nt+1,2:(Nr-1)*Nt+1) = G;
b = zeros(1,(Nr-1)*Nt+1);
b(1,1) = 1 - 4 * lambda1_r;
for i = 0:Nt-1;
    b(1,(Nr-1)*i+2) = 4*lambda1_r/Nt;
end;
H(1,1:(Nr-1)*Nt+1) = b;

r = 1/Nr:1/Nr:(Nr-1)/Nr;
fi = 1 - 2*lambda1_t./(r.*r);
ti = lambda1_t./(r.*r);
c = diag(fi);
d = diag(ti);
K = kron(I,c) + kron(J,d);
L = zeros((Nr-1)*Nt+1); L(2:(Nr-1)*Nt+1,2:(Nr-1)*Nt+1) = K;
g = zeros(1,(Nr-1)*Nt+1);
g(1,1) = 1 - 4 * lambda1_r;
for j = 0:Nt-1;
    g(1,(Nr-1)*j+2) = 4*lambda1_r/Nt;
end;
L(1,1:(Nr-1)*Nt+1) = g;

M = inv(H)*L;
N = F * M;
opts.disp = 0;
abs(eigs(N,1,'lm',opts))

eigs(N,2,'lm',opts)