// Submitted by Dr. Triranjita Srivastava. Assistant Professor, Physics Dept., Kalindi College, University of Delhi

// Aim: To prove the orthogonality of Hermitian differential Operator

// Two Hermitian Differential Operators (-id/dx) and (-d^2/dx^2) are taken as an example

// Finite Difference Method is used to formulate the matrices coresponding to the considered Differential Operator

// This method takes the value of eigenfunction equal to 0 at the initial (x0) and final boundary (xn), and x is defined as x=x0+i*h (h is step size and i is integer)
// N is the number of interval, N should be taken as odd and such that step size is small enough for high accuracy

// Using Central Differences, tridiagonal matrix is obtained for (-d/dx) which has 0 as diagonal element and -1 as upper adjacent diagonal and 1 as lower adjacent diagonal elements; 
 
//Similarly, using Central Differences, tridiagonal matrix is obtained for (-d^2/dx^2) which has 2 and -1 as upper and lower adjacent diagonal elements



clear 
clc
//********************************************************************
// Boundary over which the function is to be solved
//********************************************************************
x0=input("Enter the first boundary, x0 = ");
xn=input("Enter the second boundary, xn = ");
N = input("Enter number of intervals, N =  ");
h = (xn-x0)/N;    //step size
s=input("Enter choice 1 or 2 for differential operator -id/dx or -d^2/dx^2, respectively:  ")

select s
case 1
    //********************************************************************
    // Defining D1 Matrix corresponding to differential operator -id/dx
    //********************************************************************
    D1=zeros(N-1,N-1);

    for i=1:(N-1)
        x1(1,i)=x0+i*h;
        D1(i,i)=0;
        if i<(N-1)
            D1(i,i+1)=-%i;
            D1(i+1,i)=%i;
        end
    end
    Final_D1=D1/2*h;
    
    //*************************************************************************
    // Finding eigenvalue and eigenvector of differential operator -id/dx
    //*************************************************************************
    [eigenvector,eigenvalue] = spec(Final_D1);

case 2
    //********************************************************************
    // Defining D2 Matrix corresponding to differential operator -d^2/dx^2
    //********************************************************************

    D2=zeros(N-1,N-1);

    for i=1:(N-1)
        x1(1,i)=x0+i*h;
        D2(i,i)=2;
        if i<(N-1)
            D2(i,i+1)=-1;
            D2(i+1,i)=-1;
        end
    end
    Final_D2=D2/h^2;

    //*************************************************************************
    // Finding eigenvalue and eigenvector of differential operator -d^2/dx^2
    //*************************************************************************
    [eigenvector,eigenvalue] = spec(Final_D2);
end

    //*************************************************************************
    // Ploting of first three Eigenvector of differential operator -d^2/dx^2
    //*************************************************************************
    
    x=[x0,x1,xn];

    if s==1 then
        title('3 Lowest Order Eigenvectors of -id/dx','fontsize',4);
    else 
        title('3 Lowest Order Eigenvectors of -d^2/dx^2','fontsize',4);
    end
    
    for k = 1:3
        subplot(3,1,k)
        ylabel('A (m)','fontsize',4)
        a=get("current_axes");//get the handle of the newly created axes
        a.font_size=2
        t=get("hdl")   //get the handle of the newly created object
        t.font_size=2;  
        E_vector = [0;eigenvector(:,k);0];
        plot(x,E_vector','linewidth',2);
    end
            xlabel('x-coordinate (m)','fontsize',4)
    //*************************************************************************
    // Orthogonality Check of eigenvector of differential operator 
    //*************************************************************************
    for i=1:3
        for j=1:3
            P(i,j) = clean(sum((eigenvector(:,i).*conj(eigenvector(:,j)))));
            if i~=j & P(i,j) ~=0
                disp(" The Eigenvectors of the considered Hermitian Differential Operator are Not Orthogonal")
                abort;
            end
        end
    end
    disp(" The Eigenvectors of the considered Hermitian Differential Operator are Orthogonal")

