qm1d_fast(pot_filename,NPTS,NSTM,L)

% This routine do:
% 1. Read potential from input file.
% 2. Do spline interpolation for NPTS points.
% 3. Discretize schrodinger equation over the specified set of points:
%     v(x)=2m*V(x)/hbar^2; ee=2m*E/hbar^2.
%     [    d^2       2m        ]           2m*E
%     [ - -----  + ------ V(x) ] phi(x) = ------  phi(x)
%     [   d x^2    hbar^2      ]          hbar^2
%
% 4. Find eigen values and eigen vectors.
% 5. Plot the found eigen values and eigen vectors.
%
% Usage:
% [ee,ev] = qm1d_fast(pot_filename,NPTS,NSTM,L)
% pot_file - string that specifies filename with the potential
% NPTS - number of points for discretization of schrodinger equation
% NSTM - number of eigen values and eigen vector to find
% L - the length of the interval, starting from zero.
%
% Example: qm1d_fast('pot.dat',10000,5,10);
%
% See also: http://iffwww.iff.kfa-juelich.de/~ekoch/DFT/qm1d.html
function [ee,ev] = qm1d_fast(pot_filename,NPTS,NSTM,L)
pot_dat=load(pot_filename);
j=1:NPTS; % indexes for main diagonal
h=L/(NPTS-1); % space step

x=j*h;
V=spline(pot_dat(:,1),pot_dat(:,2),x);

main_diag=2/h^2+V(j);
sub_diag=-1/h^2*ones(1,NPTS-1);

[ee,ev] = trideigs(main_diag,sub_diag,'I',1,NSTM);

% average spacing of energy levels (for adjusting scale of ev)
de0=(ee(NSTM)-ee(1))/(NSTM-1);

% plot potential
plot(j,V(j),'r'); hold on;

% plot eign vectors
de=0.15*sqrt(NPTS)*de0;
for n=1:NSTM
    plot(j,ee(n)+de*ev(:,n)); hold on;
    plot(j,ones(length(j),1)*ee(n)); hold on;
end

% set up plotting options
xlim([0 NPTS]); ylim([min(V) ee(NSTM)+de0]);