function ffsong()  
    close all;
    dh=.1;% dh= buoc chia
    ro=1.02;% ro = mat do
    f=5500;% f =tan so
    D=100;% D= do sau
    %c=1500;la toc truyen song tai day
    N=ceil(D/dh);% so nut 
    clf;
    c=linspace(1500,1450,N)% tao mot vecto co gia tri ngau nhien
    z=dh:dh:D;% do sau, xet trong he toa do de cac
    
    % giai voi dieu kien bien tren  1
    A=make_model(dh,c,f,N); % tinh ma tran A
    [V,K]=SolveEigsA(A); % tinh tri rieng 
  
    % ve hinh 
    figure(1);
    subplot(1,4,1);% tao cac truc trong mot phan cuar so do hoa so hang va so cot duoc chia sau do chon vung ve
    plot(V(:,1:5),z);% ve do thi tuyen tinh trong khong gian 2 chieu ve gia tri x theo gia tri y
    xlabel ('v'); % chu tich van toc
    ylabel ('do sau');% chu thich do sau
    subplot(1,4,2);
    %     dua vec to V ve gia tri duong
    for m=1:5
       V(:,m)=V(:,m).*(V(:,m)./(abs(V(:,m))))
    end;
    plot(V(:,1),z,'r');
    xlabel('V mode 1');
    ylabel('do sau');
    hold on;
    [K(:)]
    % tinh truong song
    mode =5;% tinh 5 mode lon nhat cua truong song
    r = 100:1:1000;
    r=r*1e3;
    zSource = ceil(D/dh/(5/3));
    Pzr=Pressure(K(1:mode),V(:,1:mode),r,zSource);
    [zz rr]=meshgrid(z,r);
    figure(2);
    surface(rr,zz,abs(Pzr'));
    shading interp
    colorbar;
    
function A = make_model (adh,c,f,N)
    g=ones(1,N-1);
    omega=2*pi()*f;
    d=-2+adh.^2.*omega^2./c.^2.*ones(1,N);
    A=diag(d) + diag(g,1) + diag(g,-1);
       % diem o day
    A(end,end-1)=0
    A(end,end-2)=0
    
function [V,K]=SolveEigA(A)
    [V,K]=eig(A);
    K=diag(K);
    [K, index]=sort(K);
    index=flipud(index);
    K=flipud(K);
    K=sqrt(K);
    V=fliplr(V);
    
function [V,K]=SolveEigsA(A)
    [V,K]=eigs(sparse(A),50,'lm');
    K=diag(K);
    K=sqrt(K);
%    
function Pzr=Pressure(K,V,r,zSource);
    Pzk = V*diag(V(zSource,:));
    Pkr = K(:)*reshape(r,1,length(r));
    Pkr = exp(i*Pkr)./sqrt(Pkr);
    Pzr = Pzk*Pkr *i*exp(i*pi()/4)/sqrt(8*pi);
   
        