Matlab version of qm1d - 1d Schrödinger equation solver

Find eigen values and eigen vectors of Schrödinger equation and plot them.

1d Schrödinger equation solver.
This is Matlab version of qm1d program written in Fortran.
I found it here:
http://iffwww.iff.kfa-juelich.de/~ekoch/DFT/qm1d.html
Note that this program uses lapack dstevx routine, that has been wrapped by me. See my previous submit of trideigs.c.
Potential v(x) is specified with function handle to the routine parameter.

Before to run this routine install the MEX file of trideigs. See Readme-install.txt in:
http://www.mathworks.com/matlabcentral/fileexchange/24068-tridiegs/all_files

% This routine do:
% 1. Discretize Schrödinger 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
%
% 2. Find eigen values and eigen vectors.
% 3. Plot the found eigen values and eigen vectors.
%
% Usage:
% [ee,ev] = qm1d_fast(NPTS,NSTM,a,b,f_pot_handle)
% NPTS - number of points for discretization of schrodinger equation
% NSTM - number of eigen values and eigen vector to find
% a - the start point of te interval for x
% b - the end point of the interval for y
% f_pot_handle - handle to function which defines the potential
%
% Examples:
%
% 1. Harmonic oscillator - symetric at x=L/2
% NPTS=1000;
% NSTM=5;
% L=10d0;
% f_pot = @(x) {4d0*(x-L/2).^2};
% qm1d_fast(NPTS,NSTM,0,L,f_pot);
%
% 2. Harmonic oscillator - symetric at x=0
% NPTS=1000;
% NSTM=5;
% L=10d0;
% f_pot = @(x) {4d0*x.^2};
% qm1d_fast(NPTS,NSTM,-L/2,L/2,f_pot);
%
% 3. Barrier - defined with heaviside function
% NPTS=1000;
% NSTM=5;
% L=100d0;
% bb=3d0;
% aa=1d0;
% f_pot = @(x) {heaviside(x-aa)-heaviside(x-bb)};
% qm1d_fast(NPTS,NSTM,0,L,f_pot);
%
% If you don't have heaviside defined in your Matlab version define this function as:
% function [y]=heaviside(x)
% y=zeros(1,length(x))
% ipos = find(x>0);
% y(ipos)=1;
%
% 4.Barrier
% NPTS=1000;
% NSTM=5;
% L=100d0;
% bb=3d0;
% aa=1d0;
% pot_dat = [zeros(1,int16(aa/L*NPTS)), ones(1,int16((bb-aa)/L*NPTS)), zeros(1,int16((L-bb)/L*NPTS))];
% f_pot = @(x) {pot_dat(int16(x/L*(NPTS-1)))};
% qm1d_fast(NPTS,NSTM,0,L,f_pot);
%
% 5. Harmonic oscillator - with potential defined in file
% Define the potential in file:
% NPTS=1000;
% L=10d0;
% dx=L/(NPTS-1);
% pot_dat=[[0:NPTS-1]*dx; 4d0*([0:NPTS-1]-((NPTS-1)/2)).^2*dx^2]';
% save('pot.dat','-ascii','pot_dat');
% clear
%
% pot_dat=load('pot.dat');
% f_pot = @(x) {spline(pot_dat(:,1),pot_dat(:,2),x)};
% %It is not obligatory the NPTS to be the same as the number of points for the potential
% NPTS=10000;
% NSTM=5;
% L=10d0;
% qm1d_fast(NPTS,NSTM,0,L,f_pot);