function [evects,evals,Heff] = schrodinger(m,w,e)
% Jacob Biamonte (jacob.biamonte@wolfson.ox.ac.uk) 
% V 0.1
tic;
doplot = 1;
hbar = 1; 

range = 1.5; % range over which wave function is calculated --- x axis 
N = 500; %400 % For small tunnelingrates this must be increased to get reasonable answers
Neigs = 3; % Number of eigenfcns to calculate

opts.disp=0; opts.tol=1e-10; opts.maxit=10000; opts.issym=1; opts.isreal=1;

minx = - range; 
maxx = range; %
xinc = (maxx - minx) / N; % step size
x = (minx+xinc):xinc:maxx;

% define operator: 
D = sparse(1:N,1:N,-2*ones(1,N),N,N);
E = sparse(2:N,1:N-1,ones(1,N-1),N,N);
dx = 2*xinc;
S = (-1/(dx^2))*sparse(E + D + E');
% KinMat is the kinetic term in the Hamiltonian
KinMat = sparse((2*m)^(-1)*S);
% After it is calculated we wipe unnecessary variables to save memory.
clear S; clear D; clear E; 
% The N-element vector U stores the diagonal potential terms.

U = (1/2)*m*w^2*(x).^2 + e*(x);
% Energy shift to keep it positive for convenience
% U = U - min(U);

H = sparse(KinMat + sparse(1:N,1:N,U,N,N));
[evects,evals] = eigs(H,Neigs,-100,opts);
evals = diag(evals);
[evals,index] = sort(evals);
evects = (1/sqrt(xinc))*evects(:,index);

if(doplot)
    clf;
   hold on;
   % use the following to set limits of plot 
%     ylim([-2 5]);
%     xlim([minx maxx]);
    plot(x,U); 
    for cnt=1:length(evals)
       scaleevect = (evals(3) - evals(1));
       mycolor = 'k';
       plot([min(x) max(x)],[evals(cnt) evals(cnt)],[mycolor ':'])
       plot(x,((0.5*evects(:,cnt))*scaleevect + evals(cnt)),mycolor, 'lineWidth', 1)
    end
    xlabel('Position ($x$)', 'FontSize', 12, 'Interpreter', 'latex')
    ylabel('Energy', 'FontSize', 12, 'Interpreter', 'latex')
    title('Energy levels - Problem 4 sheet 2', 'FontSize', 12, 'Interpreter', 'latex')
    hold off; 
end; 











