a = 1e-3*[1 1/2 2];
ell = [1 1 2];
dx = .025;
N = ell/dx;
A3 = 2*pi*a(3)*dx;
As = 4*pi*1e-6;
rho = A3/As;
%N = [40 40 30]/10;
%dx = 1/N(1);
R2 = 0.034;
gL = 0.3;
lam = a/(2*R2*gL)/dx^2;   % lammda^2
r = a/a(3);
Hd = [2*lam(1)*ones(1,N(1)) 2*lam(2)*ones(1,N(2)) 2*lam(3)*ones(1,N(3)+1)];
Hd(1) = lam(1);
Hd(N(1)+1) = lam(2);
Hd(N(1)+N(2)+1) =  lam*r';
Hd(end) = rho*lam(3);

Hu = [-lam(1)*ones(1,N(1)-1) 0 -lam(2)*ones(1,N(2)) -lam(3)*ones(1,N(3))];
Hl = [-lam(1)*ones(1,N(1)-1) 0 -lam(2)*ones(1,N(2)-1) -r(2)*lam(2) -lam(3)*ones(1,N(3))];
Hl(end) = rho*Hl(end);

H = diag(Hd) + diag(Hu,1) + diag(Hl,-1);

H(N(1)+N(2)+1,N(1)) = -r(1)*lam(1);
H(N(1),N(1)+N(2)+1) = -lam(1);

Dd = [a(1)*ones(1,N(1)) a(2)*ones(1,N(2)) a(3)*ones(1,N(3)) a(3)/rho];
D = diag(Dd);
sD = diag((Dd).^(1/2));
sDi = diag((Dd).^(-1/2));
S = sD*H*sDi;

%norm(H'*H-H*H')   % H is not normal

%return

[Q,Z] = eig(S);

z = diag(Z);

[sz,si] = sort(z);

Qs = Q(:,si);

x3 = 0:dx:ell(3);
x1 = ell(3):dx:ell(3)+ell(1);
x2 = ell(3):dx:ell(3)+ell(2);

for k=1:9,
    q = sDi*Qs(:,k+1);
    q1 = fliplr(q(1:N(1))');
    q2 = fliplr(q(N(1)+1:N(1)+N(2))');
    q3 = fliplr(q(N(1)+N(2)+1:end)');
    subplot(3,3,k)
    stairs(x3,q3,'k')
    hold on 
    stairs(x2,[q3(end) q2],'b')
    stairs(x1,[q3(end) q1],'r')
    hold off
    axis([0 3 -10 10])
end