% Usage: [E,psi,x,Vs] = schrodinger_fd(N, 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 second-order finite-difference scheme with N points.
%
% 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 potential V(x) evaluated at x(:).

function [E,psi,x,Vs] = schrodinger_fd(N, V, m)
  x = linspace(-1,1,N+1); x = x(1:end-1);
  dx = x(2) - x(1);
  o = ones(1,N);
  A = -2 * diag(o) + diag(o(1:end-1),1) + diag(o(2:end),-1);
  A(1,N) = 1;
  A(N,1) = 1;
  Vs = V(x);
  A = - A / dx^2 + diag(Vs);
  B = diag(o);
  [psis,Es] = eig(A,B);
  E = diag(Es);
  [E,Ei] = sort(E);
  E = E(1:m);
  psi = psis(:,Ei(1:m));