function FGM
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%5
%  FGM
%       Uses Fourier Grid Method to calculate eigenstates and eigenvalues of a 
%       Hamiltonian of the form H=T+V where T, V are Kinetic and potential
%       energies respectively.
%              V_nm =  V(x_n)delta_nm
%       (d/dx)^2_nm = -k^2/3;                            for n=m
%                     -2*k^2*(-1)^(n-m)/(pi*(n-m))^2     for n!=m
%        where k=pi/dx and T=-(hbar^2/2m)(d/dx)^2;
%
%      The eigenvalues/functions are found simply by diagonalizing the matrix
%      H_nm=K_nm+V_nm 
%
% 
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

  N=2^8;          % Number of grid points
  xmax=7;         % x coordinate range
  omega=1;        % harmonic level spacing
  hbar =1;        % hbar
  M=1;            % particle mass

  dx=2*xmax/(N-1);
  x=[0:N-1]*dx-xmax;

  k=pi/dx;

  v=0.5*M*(omega*x).^2;  % harmonic potential
  
  V=diag(v);

  T=-eye(N)*k^2/3;       % kinetic energy

 
  for n=1:N
     for m=1:n-1       
        T(n,m)=-2*k^2*(-1)^(n-m)/(pi*(n-m))^2;
        T(m,n)=T(n,m);
     end;
  end;

  T=-(0.5*hbar^2/M)*T;
%keyboard;
  [U,D]=eig(T+V);

  Nplt=10;
  subplot(211);plot(x,v/8,'r',x,U(:,1:Nplt)+.25*ones(N,Nplt)*D(1:Nplt,1:Nplt));
      axis([-7,7,0,3]);
     title('HO Eigenfunctions')
     xlabel('x')
     
  subplot(212);bar(diag(D(1:Nplt,1:Nplt)),.5);
     title('HO Energy spectrum');
     xlabel('n'); ylabel('E_n');

  return;
