num_of_states=3;

n=101;
L=1.0;
pote=10;
xmax=5.0;
xmin=-1*xmax;
dx=(xmax-xmin)/(n-1);
for in=1:n
	 x(in)=xmin+dx*in;
if (abs(x(in))<L)
    pot(in)=0.0;
    else
      pot(in)=pote;
end

end


%Form the tri-diagonal matrix: assume that derivatives vahish at the edge
H=zeros(n,n);
for in=1:n
	 if (in~=1 && in~=n)
	   H(in,in)=pot(in)-0.5*(-2)/dx^2;
H(in,in+1)=-0.5/dx^2;
H(in,in-1)=-0.5/dx^2;
 else
   H(in,in)=pot(in);
end

end


%Solve n all eigenstates: probably only few lowest one reliable
%sparse solver would be better for this really, but in 
%1D we don't really care

[v,d]=eig(H);
for in=1:n
eval(in)=d(in,in);
end

%The solver might leave the solutions in a messed up order
%sort them in ascending order
[eval,iv]=sort(eval);

%Pick the corresponding eigenvectors: ground state...num_of_states
for in=1:num_of_states
mode(:,in)=v(:,iv(in));
end

%Normalize the modes (cannot trust computers to always do this by default)
for in=1:num_of_states
norm=sum(abs(mode(:,in)).^2)*dx;
mode(:,in)=mode(:,in)/sqrt(norm);
end

energies=eval(1:num_of_states)

fignum=1;
figure(fignum);clf;
fs=16;
e1=ones(size(pot))*eval(1);
e2=ones(size(pot))*eval(2);
e3=ones(size(pot))*eval(3);

subplot(2,1,1);plot(x,pot,'k','LineWidth',2);
hold on;
plot(x,e1,'b--',x,e2,'r--',x,e3,'c--','LineWidth',2);

grid on;
xlabel('x','FontSize',fs);
xlabel('E','FontSize',fs);
a=gca;set(a,'FontSize',fs);
title('Finite potential well: V=10','FontSize',fs);

axis([min(x) max(x) 0 1.1*pote]);

%for in=1:num_of_states
subplot(2,1,2);plot(x,real(mode(:,1)),'k','LineWidth',2);
grid on;hold on;
plot(x,real(mode(:,2)),'b','LineWidth',2);
plot(x,real(mode(:,3)),'r','LineWidth',2);
xlabel('x','FontSize',fs);
ylabel('Re[\psi_n(x)]','FontSize',fs);
%end
a=gca;set(a,'FontSize',fs);
hold off;
axis([min(x) max(x) -1 1]);