% Usage: [E,psi,x,Vs] = schrodinger_galerkin(N, rho, V, m)
%
% Return the m lowest eigenvalues E and eigenfunctions psi of the Schrodinger
% operator for a potential function V(x), on [-1,1] with periodic boundaries.
%
% Uses a linear-element Galerkin scheme with N points whose spacing
% is proportional to the function rho(x).  (See makegrid.m).
%
% E(i) is an array of eigenvalues, and the columns psi(:,i) are the
% eigenfunctions corresponding to E(i) at the points x(:).  Vs(:)
% is an array giving the average potential V(x) at the x(:) elements.

function [E,psi,x,Vs] = schrodinger_galerkin(N, rho, V, m)
  x = makegrid(N, rho);       % x_1 .. x_N
  xm = [ x(N)-2, x(1:N-1) ];  % x_0 .. x_{N-1}
  xp = [x(2:N), x(1)+2];      % x_2 .. x_{N+1}
  dx = xp - x(1:N);
  dxm = x(1:N) - xm;

  B = diag(dxm+dx)/3 + diag(dx(1:N-1),1)/6 + diag(dxm(2:N),-1)/6;
  B(1,N) = dxm(1)/6;
  B(N,1) = dx(N)/6;

  % check: B should be Hermitian
  if max(max(abs(B - B'))) > 1e-8
    error('B is not Hermitian!');
  end

  K = diag(1./dxm + 1./dx) - diag(1./dx(1:N-1),1) - diag(1./dxm(2:N),-1);
  K(1,N) = -1/dxm(1);
  K(N,1) = -1/dx(N);

  Vm = zeros(N,N);
  for n = 1:N
    P = quadl(@(y) V(y + xm(n)) .* y.^2, 0, dxm(n), 1e-8) / dxm(n)^2;
    M = quadl(@(y) V(y + x(n)) .* (dx(n) - y).^2, 0, dx(n), 1e-8) / dx(n)^2;
    Vm(n, n) = P + M;
  end
  for n = 1:N-1
    Vm(n,n+1) = quadl(@(y) V(y + x(n)) .* y .* (dx(n)-y), 0, dx(n), 1e-8) / dx(n)^2;
    Vm(n+1,n) = Vm(n,n+1);
  end
  Vm(N,1) = quadl(@(y) V(y + x(N)) .* y .* (dx(N)-y), 0, dx(N), 1e-8) /dx(N)^2;
  Vm(1,N) = Vm(N,1);

  A = K + Vm;
  
  % check: A should be Hermitian
  if max(max(abs(A - A'))) > 1e-8
    error('A is not Hermitian!');
  end

  [psis,Es] = eig(A,B);
  E = diag(Es);
  [E,Ei] = sort(E);
  E = E(1:m);
  psi = psis(:,Ei(1:m));

  Vs = diag(Vm) ./ diag(B); % potential values on grid