function proj632
global L T S U mul mur mu0 muL
L = 2; T=1; S=1;U=1;mul=2;mur = 10; mu0=8; muL= 1.2;
nwave = 50; nwavefac = 0.05;
%% -------------------------------------------------------------
%lambdaK = heleShawEigenSys(nwave,nwavefac);
%figure(4);
%plot((1:nwave)/nwavefac,max(lambdaK,[],1),'r*-');
% 
%% -------------------------------------------------------------
% figure(8)
% Lv = 1:0.5:10; 
% sigmax = zeros(length(Lv),1);
% ii=0;
% for muL = [2 4 5.9 7.9]
%     ii=ii+1;
%     for kk=1:length(Lv)
%         L=Lv(kk);
%         lambdaK = heleShawEigenSys(nwave,nwavefac);
%         sigmax(kk) = max(lambdaK(:));
%     end
%     subplot(2,2,ii)
%     plot(Lv,sigmax,'r*-');
%     xlabel('L');ylabel('\sigma_{max}');
%     title(['\mu(-L)=',num2str(muL)]);
% end
%% -------------------------------------------------------------
figure(9)
clear L T S U mul mur mu0 muL
T=1; S=1;U=1;mul=2;mur = 10; mu0=8; 
muLv = mul:0.83:mu0; 
Lv = [2 4 6 8]
sigmax = zeros(length(muLv),1);
ii=0;
for L = Lv;
    ii=ii+1;
    for kk=1:length(muLv)
        muL=muLv(kk)
        lambdaK = heleShawEigenSys(nwave,nwavefac);
        sigmax(kk) = max(lambdaK(:));
    end
    subplot(2,2,ii)
    plot(muLv,sigmax,'r*-');
    xlabel('\mu(-L)');ylabel('\sigma_{max}');
    title(['L = ',num2str(Lv(ii))]);
end


disp('End...');


function lambdaK = heleShawEigenSys(nwave,nwavefac)
global L T S U mul mur mu0 muL

M = 50; % # of points, #-2 interior points
d = L/(M-1);
A = zeros(M,M); B = A;

lambdaK=zeros(M,nwave);
for kk = 1:nwave
    k=kk/nwavefac;
    A(1,1) = 1/get_p(k)/d - get_q(k)/get_p(k);
    A(1,2) = -1/d/get_p(k);  % left boundary x = -L
    A(M,M) = -1/d/get_r(k) - get_s(k)/get_r(k);
    A(M,M-1) = 1/get_r(k)/d; % right boundary x = 0
    B(1,1) = 1; 
    B(M,M) = 1;
    
    for ii = 2:M-1
        x = -(ii-1)*d;
        B(ii,ii) = k^2*U*(get_mu(x+d/2)-get_mu(x-d/2))/d;        
        A(ii,ii-1) = -get_mu(x+d/2)/d^2; % grid pt: i-1/2
        A(ii,ii)   = (get_mu(x-d/2)+get_mu(x+d/2))/d^2 + get_mu(x)*k^2;
        A(ii,ii+1) = -get_mu(x-d/2)/d^2; % i+1/2
    end
    
    [f, lamda] = eig(inv(B)*A);
    lambdaK(:,kk) = 1./diag(lamda);
    tmp = 1./diag(lamda);
    %plot(tmp(tmp>
end


function y=get_mu(x)
global L T S U mul mur mu0 muL
y = (mu0-muL)/L*x+mu0;


function p=get_p(k)
global L T S U mul mur mu0 muL
p = ((mur-mu0)*U*k^2-T*k^4)/mu0;

function q=get_q(k)
global L T S U mul mur mu0 muL
q = -mur*k/mu0;

function r=get_r(k)
global L T S U mul mur mu0 muL
r = ((mul-muL)*U*k^2+S*k^4)/muL;

function s=get_s(k)
global L T S U mul mur mu0 muL
s = mul*k/muL;