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


\begin{center}
\Large \bf
Numerieke methoden in de natuurwetenschappen\\
\large 
\end{center}
\vspace{0.4cm}
\noindent 
{\bf Ex.~3a: Bound states in a finite square well, solution by use of NAG 
library}\\

Consider the problem of calculating the bound states in a potential well of 
finite height you have already solved. Solve it again by calling one of the 
routines of the NAG library for the root searching.
The NAG library is a collection of  ready-to-use subroutines which can be 
called from a f77, f90 or C program. It is useful to learn using them. 
The first step is to identify the appropriate 
subroutine, the second to understand which is the required input and 
where the output is given. A list of the avalaible subroutines can be found by 
invoking nagdtext or by looking at 
the manuals on paper which can be found in the computer room of mathematics 
on the second floor.  A description of how to use them can be found at:

http://www.sci.kun.nl/cncz/software/local/html/nag.html.

In short, you have to give the path by typing:

setenv LD\_LIBRARY\_PATH /usr/lib:/usr/ccs/lib:/usr/local/lib

and compile with the options:


f90 myprog.f90~~  -lnag -lF77

or (hopefully, let me know if it works)

CC  myprog.c 
%-L/opt/SUNWspro/lib -R/opt/SUNWspro/lib 
-lnag -lF77 -lM77 -lsunmath -lm -lc


The -lF77 option is needed because the NAG library is written in F77.

Find a subroutine which performs a root searching and insert a call to it in 
your program. Check the results against your own.
\end{document}




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 The repetitive nature of the tangent suggests that there might be 
several solutions, one per each cycle of the function. Establish the 
intervals in which roots exist. 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, Newton and secant
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}
Do not forget to label the axis 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.
%calculate the wavefunction corresponding to the first bound 
%state for two values of $V_0$ and comment the results.
%\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}
\footnote{ More details on this problem and on computational strategy can be 
found in: P.L. De Vries {\it A first course in computational Physics},
J. Wiley \& sons, N.Y. 1994}
\end{document}


