
function [A1,d]=eig_DVR(w,N)


% orthogonality relation:
% int_a^b w(x) f_n(x) f_m(x) dx = h_n \delta_{nm}
% general recursion relation:
% a_{1n} f_{n+1}(x) = (a_{2n}+a_{3n} x) f_n(x) - a_{4n} f_{n-1}(x) %rewriting the recursion relation:
% xf_n(x) = \frac{a_{1n} f_{n+1}(x) - a_{2n} f_n(x) + a_{4n} f_{n-1}(x)}{a_{3n}}
% therefore:
% \int_a^b f_{n+1}(x) x f_n(x) w(x) dx = \int_a^b \frac{a_{1n}}{a_{3n}}f^2_{n+1}(x)w(x)dx % =\frac{a_{1n}}{a_{3n}} h_{n+1}
% \int_a^b f_{n}(x) x f_n(x) w(x) dx % =- \frac{a_{2n}}{a_{3n}} h_{n}
% \int_a^b f_{n-1}(x) x f_n(x) w(x) dx % =\frac{a_{4n}}{a_{3n}} h_{n-1}
%
% a_{1n} a_{2n} a_{3n} a_{4n} h_n
% 0.Chebyshev(1st kind) t_n(x) 
% 1.Chebyshev(2nd kind) u_n(x) 1 0 2 1 pi/2
% 2. Legendre p_n(x) n+1 0 2n+1 n 2/(2n+1)
% 3.Hermite h_n(x) 1 0 2 2n pi^{1/2} 2^n n!
% 4. Laguerre l_n(x) n+1 2n+1 -1 n 1
% 5. Gen. Laguerre l_n^(a)(x) n+1 2n+a+1 -1 n+a (n+a)!/n!

%N = 12 %number of basis functions   User input
%N=50;  User input
%w=1;   User input
ip=3;
oprtr= '(w^2*X.^2)/2';
hbar=1;
m=1;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

for indx=1:N          
    n=indx-1;
    a1(indx)=1;            
    a2(indx)=0;
    a3(indx)=2;
    a4(indx)=2*n;
    h(indx)=pi^0.5*2^n*gamma(n+1);
end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% calculating the matrix X. which is x=Phi(i)*x*Phi(j) 

X(2,1) = a1(1)/a3(1)*h(2)/(h(2)*h(1))^.5;
X(1,1) = -a2(1)/a3(1);
X(N,N) = -a2(N)/a3(N);
X(N-1,N) = a4(N)/a3(N)*h(N-1)/(h(N-1)*h(N))^.5;

 for indx=3:N
     n=indx-1;
     X(indx,indx-1)= a1(n)/a3(n)*h(n+1)/(h(n+1)*h(n))^.5;
     X(indx-1,indx-1)=-a2(n)/a3(n);
     X(indx-2,indx-1)=a4(n)/a3(n)*h(n-1)/(h(n-1)*h(n))^.5;
 end


v=eval(oprtr);
[U,lamda]=eig(X);
U2=U^-1;
lamda=diag(lamda);
V=(w^2*lamda.^2)/2;
V=diag(V);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% T is a matrix for the 2nd derivative of the orthogonal polinomials. 
% T= Phi(i)*(d/dx)^2*Phi(j).  

T=0;               % matrix that represant the (d/dx)^2.
                   %  ( kinetic energy)

T(0+1,0+1) =  -0-0.5;
T(1+1,1+1) = -1-0.5;
T(N-1+1,N-1+1) = -(N-1)-0.5;
T(N-2+1,N-2+1) = -(N-2)-0.5;

T(2+1,0+1) = sqrt(2)/2;
T(3+1,1+1) = sqrt(6)/2;
T(N-3+1,N-1+1) = sqrt((N-1)*(N-1-1))/2;
T(N-4+1,N-2+1) = sqrt((N-2)*(N-3))/2;


    for indx=4:N-1       
        n=indx-2;
        T(indx-3,indx-1)=sqrt(n*(n-1))/2;
        T(indx-1,indx-1)=-n-0.5;
        T(indx+1,indx-1)=sqrt((n+1)*(n+2))/2;
    end;


T = (-hbar^2/(2*m))*(U2*T*U);

H = T+V;
[A,B] = eig(H);
[d,K]=sort(diag(B));
A1=A(:,K);
return;
