% Spectral discretization of rho^(x)^2 u_{tt}(x,t) = u_{xx}(x,t) with 
% u_x(0,t) = 0 and u_x(1,t) + u_t(1,t) = 0 (boundary damping). 
% ME October/November 2009

 N = 256;

 [D,x] = cheb01(N);      % from Trefethen, Spectral Methods in MATLAB
 [x,w] = clencurt01(N);  % from Trefethen, Spectral Methods in MATLAB

% rhosq = 2*x.^0;                      % sample mass distribution (constant)
 rhosq = (2+sin(10*pi*x));            % sample mass distribution
% rhosq = (2+1.9*sin(10*pi*x))/4;     % sample mass distribution; use large N

 M = w*sqrt(rhosq);
 rho1 = sqrt(rhosq(1));

 asymptote = -1/(2*M)*log(abs((1+rho1)/(1-rho1)));
 fprintf('predicted asymptote = %10.7f\n', asymptote)

 figure(1), clf
 plot(x,sqrt(rhosq),'b-','linewidth',2)
 xlabel('$x$','interpreter','latex')
 ylabel('$\rho(x)$','interpreter','latex')

 L = diag(1./rhosq)*D*D;

 dNN = D(1,1);
 dNt = D(1,2:N);
 dN0 = D(1,N+1);
 dl  = D(2:N,1);
 Dhat = D(2:N,2:N);
 dr  = D(2:N,N+1);
 d0N = D(N+1,1);
 d0t = D(N+1,2:N);
 d00 = D(N+1,N+1);

 lNN = L(1,1);
 lNt = L(1,2:N);
 lN0 = L(1,N+1);
 ll  = L(2:N,1);
 Lhat = L(2:N,2:N);
 lr  = L(2:N,N+1);
 l0N = L(N+1,1);
 l0t = L(N+1,2:N);
 l00 = L(N+1,N+1);

 A = [(d0N*dN0/d00-dNN)   (d0t*dN0/d00-dNt)   zeros(1,N);
       zeros(N-1,1)   zeros(N-1)    zeros(N-1,1)  eye(N-1);
       1 zeros(1,N-1) zeros(1,N);
       ll-lr*d0N/d00  Lhat-lr*d0t/d00 zeros(N-1,N)];

 B = eye(2*N);
 B(N+1,1:N+1) = [d0N*rhosq(N+1)*w(N+1)/d00-rhosq(1)*w(1) ...
                 rhosq(N+1)*w(1)*d0t/d00-rhosq(2:N)'.*w(2:N) 0];

 [V,Lam] = eig(A,B);
 ew = diag(Lam); 


 figure(2), clf
 plot([0 0],[-10000 10000],'-','color',.7*[1 1 1]), hold on
 plot([-100 100],[0 0],'-','color',.7*[1 1 1])
 plot(asymptote*[1 1],[-10000 10000],'-','color',[1 0 0]), hold on
 plot(real(ew),imag(ew),'k.','markersize',10)
 title(sprintf('asymptote at $%10.7f$',asymptote),'interpreter','latex')

 indx = find((imag(ew)>=0)&(abs(ew)<1e8)&(real(ew)>-1000));
 ew = ew(indx);
 V = V(:,indx);
 [junk,indx] = sort(imag(ew));
 ew = ew(indx);
 V  = V(:,indx);
 ew = ew(1:floor(N/4));
 axis([-5 1 -1.25*max(imag(ew)) 1.25*max(imag(ew))])

 oplt = [];
 for k=1:min(length(ew),50)
     uN = V(1,k);
     u  = V(2:N,k);
     u0 = (-d0N*uN-d0t*u)/d00;
     U = [uN;u;u0];
     U = U/max(abs(U));
     figure(2)
     if ~isempty(oplt), delete(oplt), end
     hold on, oplt=plot(real(ew(k)),imag(ew(k)),'ro','markersize',8);
     figure(3), clf
     plot(x,real(U),'b-','linewidth',2)
     hold on
     plot(x,abs(U),'r-','linewidth',1)
     plot(x,-abs(U),'r-','linewidth',1)
     axis([0 1 -1.1 1.1])
     title(sprintf('$\\lambda=%10.7f+%10.7fi$',real(ew(k)),imag(ew(k))),...
            'interpreter','latex')
   pause
 end
 if ~isempty(oplt), delete(oplt), end