%{
  Calculation of the eigenmodes of a vibrating string.
  Numerical method: finite-difference scheme, second-order central.

  clear all
  close all
  clc

  N = 50;
  [S,V_an,kappa_an]=eigenProblem(N,{'fd2','fd4'});
%}
function [S,V_an,kappa_an]=eigenProblem(N,method)

% Analytical eigenvalues/eigenvectors
kappa_an = (1:N-1).*pi;
V_an     = sin(repmat((1:N-1),[N+1,1]).*repmat((0:N).',[1,N-1]).*(pi/N));
% Normalization
V_an     = V_an./(repmat(sqrt(sum(V_an.^2,1)),[N+1,1]));

color = {'k','g'};

for k=1:length(method)

    % Numerical eigenvalues/eigenvectors
    switch method{k}
        case{'fd2'}
            A = fd2(N-1);
        case{'fd4'}
            A = fd4(N-1);
    end

    [S(k).V,S(k).D] = eigs(A,N-1,'smallestabs');
    S(k).V = [zeros(1,N-1);S(k).V;zeros(1,N-1)];
    S(k).kappa = N*sqrt(-diag(S(k).D));

    % Reverse sign of numerical eigenvectors if necessary
    for j=1:N-1
        ep = norm(S(k).V(:,j)-V_an(:,j));
        em = norm(S(k).V(:,j)+V_an(:,j));
        if (em<ep)
            S(k).V(:,j)= -S(k).V(:,j);
        end
    end



    m = (1:N-1)./N;

    if ~exist('hf','var')
        hf=figure;
        sb1=subplot(1,2,1);
        hold on
        sb2=subplot(1,2,2);
        hold on
    end

    figure(hf);
    subplot(sb1);
    h=plot(m,kappa_an,'ro');
    h=plot(m,S(k).kappa   ,'x');
    set(h,'markeredgecolor',color{k});

    figure(hf);
    subplot(sb2);
    h=plot((0:N)./N,V_an(:,1),'r-');
    h=plot((0:N)./N,S(k).V(:,1),'k-');
    set(h,'color',color{k});

    h=plot((0:N)./N,V_an(:,round(N/4)),'r-');
    h=plot((0:N)./N,S(k).V(:,round(N/4)),'kx');
    set(h,'color',color{k});

end

figure(hf);
subplot(sb1);
box on
xl=xlabel('$m/N$','interpreter','latex');
yl=ylabel('$\widetilde{k}$','interpreter','latex');

figure(hf);
subplot(sb2);
box on
xl=xlabel('$\xi$','interpreter','latex');
yl=ylabel('$\mathbf{Y}^{(m)}$','interpreter','latex');

end

function A=fd2(dim)

one = ones(dim,1);
A = spdiags([[one(1:end-1);0],-2*one,[0;one(1:dim-1)]],[-1,0,1],dim,dim);

end

function A=fd4(dim)

one = ones(dim,1);
A = spdiags([[-one(1:end-2)./12;0;0],[one(1:end-1).*(4/3);0],-2.5*one,...
    [0;one(1:end-1).*(4/3)],[0;0;-one(1:end-2)./12]],...
    [-2:2],dim,dim);

A(1,1:2) = [-2 , 1]; A(1,3) = 0;
A(end,end-1:end) = [1 , -2]; A(end,end-2) = 0;

end