% Integrate Schrodinger's equation once for a particular potential V, given
% as a function handle, a domain width w (centered at x=0), and an assumed
% (possibly incorrect) eigenvalue E. The procedure is:
%  - Start with a V=0 boundary condition at x = -w/2.
%  - Integrate to find psi(x) on  -w/2 < x < w/2.
%  - Normalize psi, dpsi/dx so that max(psi) = 1 (nice for plotting).
%  - Evaluate deviation from V=0 boundary condition at x = w/2.
%  - Return this deviation as the error, 'err'.
% The roots of err as a function of E are the correct energy eigenvalues. 
% The returned q is an array of psi and d/dx(psi) 
% CCK 2009-Apr-05
function [err, q, x] = oneshot(E, V, w, Npoints)

if !exist('Npoints','var') 
   Npoints = 1000; % default value
end

x = linspace(-w/2, w/2, Npoints);

% boundary conditions
psi_0 = 1; % completely arbitrary
dpsi_dx_0 = psi_0 * sqrt(-E); 
   % consistent with psi --> 0 as x --> -infinity, if V=0.
q_0 = [psi_0, dpsi_dx_0];

% integrate the ODE
q = lsode(@(q,x) schrodinger(q,x,E,V), q_0, x);

% Normalize the wave function to a maximum of 1.
normalization = max( abs(q(:,1)) );
q = q/normalization;

% evaluate the error
psi_final = q(Npoints,1);
dpsi_dx_final = q(Npoints,2);
err = dpsi_dx_final + sqrt(-E) * psi_final;
   % err == 0 is equivalent to psi --> 0 as x --> infinity with V=0.
