\documentclass[12pt]{article}
\textwidth 16cm
\textheight 23cm
\topmargin 0cm
\evensidemargin 0cm
\oddsidemargin 0cm
\parindent=0pt
\parskip=10pt
\begin{document}


\begin{center}
\Large \bf
Computational Methods and Data Analysis 2003\\
\large 
\end{center}
\vspace{0.4cm}
\noindent 
{\bf Ex.~4: Bound states in a finite square well}\\

This exercise requires the implementation of algorithms for
finding the zero's of nonlinear equations and for numerical integration. 

Consider a square potential well of finite height V$_0$
and width $2a$, as in figure ~\ref{fig:put}. 


\begin{figure}[hbt]
 \centering
\setlength{\unitlength}{0.00050000in}%
%
\begin{picture}(6066,5493)(2968,-6700)
\thicklines
\put(6001,-6361){\line( 0, 1){4800}}
\put(3001,-2761){\line( 1, 0){1200}}
\put(4201,-2761){\line( 0,-1){3600}}
\put(4201,-6361){\line( 1, 0){3600}}
\put(7801,-6361){\line( 0, 1){3600}}
\put(7801,-2761){\line( 1, 0){1200}}
\put(6001,-1561){\line( 0, 1){ 75}}
\put(6001,-1486){\vector( 0, 1){ 75}}
\put(6001,-6661){\small 0}
\put(3901,-6661){\large $-a$}
\put(7801,-6661){\large $a$}
\put(6226,-1411){\large $V(x)$}
\put(6301,-2761){\large $V_{0}$}
\end{picture}
 \caption[1]{Finite potential well}
 \label{fig:put}

\end{figure}

The Schr\"odinger equation:
\[
\frac{d^2\Psi}{dx^2}-\frac{2m}{\hbar^2}\left[V(x)-E\right]\Psi=0
\]

with $V(x)$ given by:
\begin{eqnarray*}
V(x)&=&V_0 \qquad x<-a \mbox{~~~~~~~~~~~~~~~~~~~~~region $I$~~~}\\
V(x)&=&0 \qquad -a<x<a \mbox{~~~~~~~~~~~~~~~~region $II$~~}\\
V(x)&=&V_0 \qquad x>a \mbox{~~~~~~~~~~~~~~~~~~~~~~~~region $III$~}
\end{eqnarray*}

has general solutions for the bound states at energy less than $V_0$
in the three regions given by:
\begin{eqnarray*}
\Psi_{I}(x)&=&A~\exp^{\beta x}+B~\exp^{-\beta x} \qquad x<-a\\
\Psi_{II}(x)&=&C~\sin\alpha x+D~\cos\alpha x \qquad -a<x<a\\
\Psi_{III}(x)&=&E~\exp^{\beta x}+F~\exp^{-\beta x} \qquad x>a\\[1ex]
\alpha&=&\sqrt{\frac{2mE}{\hbar^2}}\\
\beta&=&\sqrt{2m(V_0-E)/\hbar^2}
\end{eqnarray*}

The boundary conditions require continuity of the wavefunction and of 
its derivative at the boundary, whence   
\\[1ex]
Even states: $B=E=0, C=0, D\neq0, A=F \qquad \alpha \tan \alpha a = \beta$\\
Odd states: $B=E=0 , C\neq0, D=0, A=-F \qquad \alpha \cot \alpha a = -\beta$\\

Therefore eigenvalues of the Schr\"odinger equation are found by solving
numerically:
\begin{center}
\begin{tabular}{ccc}
$\alpha tg(\alpha a)-\beta$=0 &for even states\\
$\alpha ctg(\alpha a)+\beta$=0 &for odd  states
\end{tabular} 
\end{center}

A more convenient form is:

\begin{center}
\begin{tabular}{ccc}
$\alpha sin(\alpha a)-\beta cos(\alpha a)$=0 &for even states\\
$\alpha cos(\alpha a)+\beta sin(\alpha a)$=0 &for odd  states
\end{tabular} 
\end{center}


Use $2a$=10~\AA, $m=1~m_e$. 
It is useful to use $\hbar^2=7.6199682 m_e$~eV~\AA$^2$   

\begin{enumerate}
\item  Bracket the range in which solutions exist for even and odd states.
This can be done by noting that $\beta/\alpha$ is always positive so that, for 
instance, solutions for the even state will occur in ranges of $\alpha a$ 
such that the tangent is positive.  
Remember that, in one dimension,  a square well
has always at least one bound state at $E\neq 0$.
\item Find all bound states with precision 10$^{-4}$~eV
by use of the bisection and Newton-Raphson
methods for several values of $V_0$ between $V_0$=0.2~eV and $V_0$=5~eV. 
\item 
Make a plot of the 
eigenvalues as a function of $V_0$, comparing also to the limiting case of an 
infinite well: 

\begin{equation}
E_n= {\hbar^2\over 2m_e}({n\pi\over 2a})^2. 
\end{equation}
Give some of the data also in a Table in order to check that the results are 
correct within the desired precision. Do not forget to label the axis of the 
figures and to make clear table and figure captions, giving the used units 
and for which values of 
the parameters  (e.g which value of $a$) the calculation is done. 

The idea is that, even without the text of the exercise, by reading the 
report one can understand and reproduce the calculation. 
Comment briefly the results.
%If you use gnuplot, you could plot the eigenvalues with linespoints so that you% can trace back for which values of $V_0$ the calculations has been performed.
\item Facultative: Make a plot of the eigenvalues as a function of 
$a$ for $V_0$=3 eV.

\item Calculate (and plot) the normalized wavefunction corresponding to the 
first two bound states for two values of $V_0$ and comment the results.
To perform the integral needed to normalize the wavefunction you need to 
specify the range of integration. This obviously requires an approximation, 
the desired range is in principle from minus to plus infinity.  However, 
once you know the energy of a given state you can estimate the decay length 
of the wavefunctions in the forbidden barrier regions in order to define 
a sensible range for integration. Describe your reasoning in doing this.
You can either perform the integrals by use of the built-in matlab functions or by constructing your own code for trapezoidal and Simpson integration. 
In defining the grid of points on which 
the wavefunction is calculated to perform the integration choose a
spacing which fits into the well. 

\item Check numerically that, as expected from quantum mechanics, 
the eigenfunctions corresponding to different eigenvalues are orthogonal. 

\item Facultative: Integrate the wavefunction of one of the states 
by means of Monte Carlo integration by writing
\begin{equation}
\int_a^b{f(x) dx}= (b-a) {1\over N} \Sigma_{i=1}^N f(x_i)
\end{equation}
where $x_i$ are $N$ random numbers uniformely distributed between $a$ and $b$.
Study the dependence on $N$.

%\item Could you find out which is the value of $V_0$ which gives the first 3
%bound states at energy $E_1$=0.26230 $E_2$=1.03633  $E_3$= 2.27006 eV?
\end{enumerate}
\end{document}


