C      ALGORITHM 810, COLLECTED ALGORITHMS FROM ACM.
C      THIS WORK PUBLISHED IN TRANSACTIONS ON MATHEMATICAL SOFTWARE,
C      VOL. 27,NO. 2,      June, 2001, P.  143--192.
#! /bin/sh
# This is a shell archive, meaning:
# 1. Remove everything above the #! /bin/sh line.
# 2. Save the resulting text in a file.
# 3. Execute the file with /bin/sh (not csh) to create the files:
#	Doc/
#	Doc/README.first
#	Doc/autoinput.txt
#	Doc/help.tex
#	Doc/intro.tex
#	Doc/readme.txt
#	Doc/xamples.tex
#	Fortran/
#	Fortran/Sp/
#	Fortran/Sp/Drivers/
#	Fortran/Sp/Drivers/data2
#	Fortran/Sp/Drivers/driver1.f
#	Fortran/Sp/Drivers/driver2.f
#	Fortran/Sp/Drivers/driver3.f
#	Fortran/Sp/Drivers/makepqw.f
#	Fortran/Sp/Drivers/res1
#	Fortran/Sp/Drivers/res2
#	Fortran/Sp/Drivers/res3
#	Fortran/Sp/Drivers/xamples.f
#	Fortran/Sp/Src/
#	Fortran/Sp/Src/src.f
# This archive created: Tue Nov 13 09:51:12 2001
export PATH; PATH=/bin:$PATH
if test ! -d 'Doc'
then
	mkdir 'Doc'
fi
cd 'Doc'
if test -f 'README.first'
then
	echo shar: will not over-write existing file "'README.first'"
else
cat << "SHAR_EOF" > 'README.first'
In order to make this code compatible with the other CALGO codes
there has been a renaming of filenames. This hasn't gone through to
the documentation; the mappings are

coupdr -> driver1
drive  -> driver2
sepdr  -> driver3
sleign2.f -> src.f

data2 needs to be renamed auto.in before running driver2 in
automatic mode.
SHAR_EOF
fi # end of overwriting check
if test -f 'autoinput.txt'
then
	echo shar: will not over-write existing file "'autoinput.txt'"
else
cat << "SHAR_EOF" > 'autoinput.txt'
		SLEIGN2: THE AUTOINPUT.TXT FILE

01 March 2001: P.B. BAILEY, W.N. EVERITT AND A. ZETTL

     The experienced user of SLEIGN2 who would like to avoid the
"question & answer" format in favor of a more direct (and faster) system
can do so by simply creating a very brief text file containing the
parameters of his particular problem. The file must be called
   
                   auto.in

and must reside in the same directory as the executable SLEIGN2 file
(possibly xamples.x or bloggs.x) .  Once such a file has been created,
the user simply types

        xamples.x  <ENTER>

as usual (or bloggs.x <ENTER>), whereupon the prompt 

   WOULD YOU LIKE AN OVERVIEW OF HELP ? (Y/N) (h?)

appears -- as usual.  But instead of replying to the question
with y, or n, or h, one simply types in the response

           a  <ENTER>

The "a" here stands for "automatic" operation of SLEIGN2;
at this point the computation of the requested eigenvalue(s) 
proceeds without further action by the user.
     The construction of the file "auto.in" consists of simply
defining a number of "KEYWORDS", each on a separate line, which
together constitute a complete set of input parameters defining
the eigenvalue problem to be solved. (The differential equation
coefficients have, presumably, been already defined in either
xamples.f, or bloggs.f, or.. .)   
          -------------------------------------------
     The KEYWORDS, all of which end in a colon and must be
followed by at least one space, are:

   a:  -- The left end-point of the interval (a,b); 
          Value is any real number;
          Default value is -infinity.
   b:  -- The right end-point of the interval (a,b);
          Value is any real number;
          Default value is +infinity.
   classa:  -- The end-point classification of a;
               Value is one of { r, wr, lcno, lco, lp, d }.
   classb:  -- The end-point classification of b:
               Value is one of { r, wr, lcno, lco, lp, d }.
   bca: -- boundary condition for the end-point a;
           Value is either d (for Dirichlet),
                           n (for Neuman),
                        or two real numbers A1, A2 . 
   bcb: -- boundary condition for the end-point b;
           Value is either d (for Dirichlet),
                           n (for Neuman),
                        or two real numbers B1, B2 . 
   bcc:  -- coupled boundary condition;
            Value is either p (for periodic),
                            s (for semi-periodic),
                        or five real numbers
                           alpha, k11, k12, k21, k22 .
   numeig:  -- Index (or range of indices) of desired eigenvalue;
               Value is an integer N1, or pair of integers N1,N2 .
   param:  -- Parameter(s) appropriate for the problem;
              Value is one or two real numbers, param1, param2.
   np:  --  Problem number;
            (Appropriate only if one of the differential equations
              in xamples.f is being used.);
            Value is an integer from 1 to 32.
   output:  -- Name of output file;
               Value is a character string, the name of the output
                 file where the results of the computation are to
                 be written;
               Default value is "auto.out" .
   end:  -- Last line of file "auto.in" ; no value set.
            
   again:  --  Terminates the input for one eigenvalue problem and
                 begins the input for another ; no value set.

Although the KEYWORDS used to define any one problem can be defined
in any order, there are a few rules to be observed:  Namely,

  Only those KEYWORDS whose values are necessary need be mentioned.
  Any KEYWORD definition remains in effect until redefined;
  To erase a previous definition of a KEYWORD, redefine it to 
     have the value "null".  
          -----------------------------------------------
     
     Perhaps a simple example would help to make clear what such
a file would look like.  One such could be:

output: bessel.rep
np: 2
param: 0.75
a: 0.0
b: 1.0
classa: lcno
classb: r
bca: 1.0,0.0
bcb: d
numeig: 2,5
end:

     This file (which must be called "auto.in", of course), would
be suitable for running Bessel's equation in xamples.f.   
Evidently the problem selected is #2 of xamples.f (Bessel's equation),
with the parameter nu = 0.75, on the interval (0.0,1.0) .  The
end-point a is asserted to be LCNO; end-point b is R; the boundary
condition at a is defined by A1 = 1.0 & A2 = 0.0; boundary condition
at b is Dirichlet; and the eigenvalues lambda(2), lambda(3), lambda(4),
lambda(5) are to be computed.
 
     If the user wishes to run several problems, one after another, 
things can become a little complicated, but can be done.  Here is
a more complicated example of a file auto.in :
 
output: Legendre 
np: 1
a: -1.0
classa: lcno
bca: 1.0,0.0
b: 1.0
classb: lcno
bcb: 1.0,0.0
numeig: 0,5
 
again:
np: 6
output: Sears-Titchmarsh
a: 1.0
classa: r
classb: lco
bcb: 1.0,0.0
numeig: -1,2 
bca: d
b: null
 
again:
numeig: 1,5
np: 12 
a: 0.0 
output: Mathieu.rep
classa: r
classb: r
b: pi
bcc: p 
param: 5.0
 
again:
np: 16
output: Jacobi.rep
a: -1.0 
b:  1.0
param: 0.5,-1.2
classa: lp
classb: lcno
bcb:  1.0,0.0
numeig: 0,2
 
again:
param: 1.2,-0.5
classa: wr
classb: lp
bca: d
bcb: null
 
end:
      ---------------------------------------------------
 
     The user must be  aware, however, that this "automatic" system
for running SLEIGN2 has no built in safety features.  
Whereas the "question & answer" format can and does catch many kinds
of simple errors on the part of the user, there are no safety
devices of any kind in place when SLEIGN2 is run in automatic mode.
If the program runs, and if an output file has been written, then
the user has at least a record of the parameters of the problem which
were used for the run.  However if the auto.in file is constructed with
any errors in it, almost anything can happen. 
   
Department of Mathematical Sciences, Northern Illinois University,
DeKalb, IL 60115-2888, USA.


SHAR_EOF
fi # end of overwriting check
if test -f 'help.tex'
then
	echo shar: will not over-write existing file "'help.tex'"
else
cat << "SHAR_EOF" > 'help.tex'
%% This document created by Scientific Word (R) Version 3.5

\documentclass[12pt]{amsart}%
\usepackage{graphicx}
\usepackage{amscd}
\usepackage{amsmath}%
\usepackage{amsfonts}%
\usepackage{amssymb}
%TCIDATA{OutputFilter=latex2.dll}
%TCIDATA{CSTFile=amsartci.cst}
%TCIDATA{Created=Tue Sep 07 11:15:04 1999}
%TCIDATA{LastRevised=Sunday, April 01, 2001 11:19:55}
%TCIDATA{<META NAME="GraphicsSave" CONTENT="32">}
%TCIDATA{<META NAME="DocumentShell" CONTENT="Journal Articles\amsnumb">}
%TCIDATA{Language=British English}
\newtheorem{theorem}{Theorem}
\theoremstyle{plain}
\newtheorem{acknowledgement}{Acknowledgement}
\newtheorem{algorithm}{Algorithm}[section]
\newtheorem{axiom}{Axiom}[section]
\newtheorem{case}{Case}[section]
\newtheorem{claim}{Claim}[section]
\newtheorem{conclusion}{Conclusion}[section]
\newtheorem{condition}{Condition}[section]
\newtheorem{conjecture}{Conjecture}[section]
\newtheorem{corollary}{Corollary}[section]
\newtheorem{criterion}{Criterion}[section]
\newtheorem{definition}{Definition}[section]
\newtheorem{example}{Example}[section]
\newtheorem{exercise}{Exercise}[section]
\newtheorem{lemma}{Lemma}[section]
\newtheorem{notation}{Notation}[section]
\newtheorem{problem}{Problem}[section]
\newtheorem{proposition}{Proposition}[section]
\newtheorem{remark}{Remark}[section]
\newtheorem{solution}{Solution}[section]
\newtheorem{summary}{Summary}[section]
\numberwithin{equation}{section}
\numberwithin{theorem}{section}
\newcommand{\thmref}[1]{Theorem~\ref{#1}}
\newcommand{\secref}[1]{\S\ref{#1}}
\newcommand{\lemref}[1]{Lemma~\ref{#1}}
\setlength{\oddsidemargin}{0.0in}
\setlength{\evensidemargin}{0.0in}
\setlength{\textwidth}{6.5in}
\setlength{\textheight}{8.5in}
\setlength{\headsep}{0.25in}
\setlength{\headheight}{0.0in}


\begin{document}
\title[HELP file]{Sleign2: the HELP file}
\author{P.B. Bailey}
\author{W.N. Everitt}
\author{A. Zettl}
\address{Department of Mathematical Sciences, Northern Illinois University, DeKalb, IL
60115-2888, USA}
\date{01 March 2001 (File: help.tex)}
\maketitle


This copy of the HELP data has been written in AMS-LaTeX in view of the
considerable amount of mathematical formulae involved in the text. The user
should note that when HELP is accessed in MAKEPQW and/or DRIVE the
corresponding data is given in Fortran notation. The best use of this HELP
data can be made by printing out a copy of this AMS-LaTeX file to have
available when the code files are in use.

HELP may be called at any point where the program halts and displays (h?), by
pressing ``h
%TCIMACRO{\TEXTsymbol{<}}%
%BeginExpansion
$<$%
%EndExpansion
ENTER%
%TCIMACRO{\TEXTsymbol{>}}%
%BeginExpansion
$>$%
%EndExpansion
''. To RETURN from HELP, press ``r
%TCIMACRO{\TEXTsymbol{<}}%
%BeginExpansion
$<$%
%EndExpansion
ENTER%
%TCIMACRO{\TEXTsymbol{>}}%
%BeginExpansion
$>$%
%EndExpansion
''. To QUIT at any program halt, press ``q
%TCIMACRO{\TEXTsymbol{<}}%
%BeginExpansion
$<$%
%EndExpansion
ENTER%
%TCIMACRO{\TEXTsymbol{>}}%
%BeginExpansion
$>$%
%EndExpansion
''.

This AMS-LaTeX file is supplied as a separate text file within the SLEIGN2
package; it can be accessed on-line in both the MAKEPQW (if used) and DRIVE files.

HELP contains information to aid the user in entering data: on the coefficient
functions $p,q,w;$ on the self-adjoint, separated and coupled, regular and
singular, boundary conditions; on the limit-circle boundary condition
functions $u,v$ at the endpoint $a$ and $U,V$ at the endpoint $b$ of the
interval $(a,b)$; on the endpoint classifications of the differential
equation; on DEFAULT entry; on eigenvalue indexes; on IFLAG information; and
on the general use of the program SLEIGN2.

\medskip

The 17 sections of HELP are:

H1: Overview of HELP. H2: File name entry. H3: The differential equation. H4:
endpoint classification. H5: DEFAULT entry. H6: Self-adjoint limit-circle
boundary conditions. H7: General self-adjoint boundary conditions. H8:
Recording the results. H9: Type and choice of interval. H10: Entry of
endpoints. H11: endpoint values of $p,q,w$. H12: Initial value problems. H13:
Indexing of eigenvalues. H14: Entry of eigenvalue index, initial guess, and
tolerance. H15: IFLAG information. H16: Plotting. H17: Indexing of eigenvalues.

HELP can be accessed at each point in MAKEPQW and DRIVE where the user is
asked for input, by pressing ``h
%TCIMACRO{\TEXTsymbol{<}}%
%BeginExpansion
$<$%
%EndExpansion
ENTER%
%TCIMACRO{\TEXTsymbol{>}}%
%BeginExpansion
$>$%
%EndExpansion
''; this places the user at the appropriate HELP section. Once in HELP, the
user can scroll the further HELP sections by repeatedly pressing ``h
%TCIMACRO{\TEXTsymbol{<}}%
%BeginExpansion
$<$%
%EndExpansion
ENTER%
%TCIMACRO{\TEXTsymbol{>}}%
%BeginExpansion
$>$%
%EndExpansion
'', or jump to a specific HELP section Hn (n =1,2,...17) by typing ``hn
%TCIMACRO{\TEXTsymbol{<}}%
%BeginExpansion
$<$%
%EndExpansion
ENTER%
%TCIMACRO{\TEXTsymbol{>}}%
%BeginExpansion
$>$%
%EndExpansion
''; to return to the place in the program from which HELP is called, press
``r
%TCIMACRO{\TEXTsymbol{<}}%
%BeginExpansion
$<$%
%EndExpansion
ENTER%
%TCIMACRO{\TEXTsymbol{>}}%
%BeginExpansion
$>$%
%EndExpansion
''.

\medskip

H2: File name entry.

MAKEPQW is used to create a Fortran file containing the coefficients $p,q,w$
defining the differential equation, and the boundary condition functions $u,v$
and $U,V$ if required. The file must be given a NEW filename which is
acceptable to your Fortran compiler. For example, it might be called bessel.f
or bessel.for depending upon your compiler.

The same naming considerations apply if the Fortran file is prepared other
than with the use of MAKEPQW.

\medskip

H3: The differential equation.

The prompt ``Input $p$ (or $q$ or $w$) ='' requests you to type in a Fortran
expression defining the function $p$, which is one of the three coefficient
functions defining the Sturm-Liouville differential equation%
\begin{equation}
-(py^{\prime})^{\prime}+qy=\lambda wy \tag{*}%
\end{equation}
to be considered on some interval $(a,b)$ of the real line. The actual
interval used in a particular problem can be chosen later, and may be either
the whole interval $(a,b)$ where the coefficient functions $p,q,w$ are
defined, or on any sub-interval $(a^{\prime},b^{\prime})$ of $(a,b)$;
$a=\infty$ and/or $b=+\infty$ are allowable choices for the endpoints.

The coefficient functions $p,q,w$ of the differential equation may be chosen
arbitrarily but must satisfy the following conditions:

(1) $p,q,w$ are real-valued throughout $(a,b)$.

(2) $p,q,w$ are piece-wise continuous and defined throughout the interior of
the interval $(a,b)$.

(3) $p$ and $w$ are strictly positive in $(a,b)$.

\noindent For better error analysis in the numerical procedures, condition (2)
above is often replaced with

(2$^{\prime}$) $p,q,w$ are four times continuously differentiable on $(a,b)$.

\noindent The behaviour of $p,q,w$ near the endpoints $a$ and $b$ is critical
to the classification of the differential equation (see H4 and H11).

\medskip

H4: endpoint classification.

The correct classification of the endpoints $a$ and $b$ is essential to the
working of the SLEIGN2 program. To classify the endpoints, it is convenient to
choose a point $c$ in $(a,b)$; \textit{i.e.} $a<c<b$. \quad\quad$\quad
\;\;$Subject to the general conditions on the coefficient functions $p,q,w$
(see H3):

(1) a is REGULAR (say R) if $-\infty<a,$ $p,q,w$ are piece-wise continuous on
$[a,c]$, and $p(a)>0,w(a)>0.$

(2) $a$ is WEAKLY REGULAR (say WR) if $-\infty<a,$ $a$ is not R and%
\[
\int_{a}^{c}\left\{  p^{-1}+\left|  q\right|  +w\right\}  <+\infty.
\]

If endpoint $a$ is neither R nor WR, then $a$ is SINGULAR; that is, either
$-\infty=a,$ or $-\infty<a$ and%
\[
\int_{a}^{c}\left\{  p^{-1}+\left|  q\right|  +w\right\}  =+\infty.
\]

(3) The SINGULAR endpoint $a$ is LIMIT-CIRCLE NON-OSCILLATORY (say LCNO) if
for some real value of the spectral parameter $\lambda$ ALL real-valued
solutions $y$ of the differential equation%
\[
-(py^{\prime})^{\prime}+qy=\lambda wy\;\text{on}\;(a,c]
\]

satisfy the conditions:%
\[
\int_{a}^{c}w\left|  y\right|  ^{2}<+\infty
\]

\noindent and $y$ has at most a finite number of zeros in $(a,c].$

(4) The SINGULAR endpoint a is LIMIT-CIRCLE OSCILLATORY (say LCO) if for some
real $\lambda$ ALL real-valued solutions of the differential equation (*)
satisfy the conditions:%
\[
\int_{a}^{c}w\left|  y\right|  ^{2}<+\infty
\]

\noindent and $y$ has an infinite number of zeros in $(a,c].$

(5) The SINGULAR endpoint a is LIMIT POINT (say LP) if for some real $\lambda$
at least one solution of the differential equation $(\ast)$ satisfies the
condition:%
\[
\int_{a}^{c}w\left|  y\right|  ^{2}=+\infty.
\]

There is a similar classification of the endpoint $b$ into one of the five
distinct cases R, WR, LCNO, LCO, LP.

Although the classification of singular endpoints invokes a real value of the
parameter $\lambda$, this classification is invariant in $\lambda$; all real
choices of $\lambda$ lead to the same classification.

In determining the classification of singular endpoints for the differential
equation $(\ast)$, it is often convenient to start with the choice $\lambda=0$
in attempting to find solutions (particularly when $q=0$ on $(a,b)$); however,
see example 7 below.

See H6 on the use of maximal domain functions to determine the classification
at singular endpoints.

EXAMPLES:

1. $-y^{\prime\prime}=\lambda y$ is R at both endpoints of $(a,b)$ when $a$
and $b$ are finite.

2. $-y^{\prime\prime}=\lambda y$ on $(-\infty,+\infty)$ is LP at both endpoints.

3. $-\left(  x^{1/2}y^{\prime}(x)\right)  ^{\prime}=\lambda x^{-1/2}y(x)$ for
all$\;x\in(0,+\infty)$ is WR at $0$ and LP at $+\infty$ (take $\lambda=0$ in
$(\ast)$). See the file xamples.f, \#10 (Weakly Regular).

4. $-((1-x^{2})y^{\prime}(x))^{\prime}=\lambda y(x)\;$for all$\;x\in(-1,+1)$
is LCNO at both ends (take $\lambda=0$ in $(\ast)$). See xamples.f, \#1 (Legendre).

5. $-y^{\prime\prime}(x)+Cx^{-2}y(x)=\lambda y(x)\;$for all$\;x\in(0,+\infty)$
is LP at $+\infty;$ however at $0$ the equation illustrates a number of
different classifications(take $\lambda=0$ in $(\ast)$) as follows:

LP $C\geq3/4$ ; LCNO for $-1/4\leq C<3/4$ (but $C\neq0$); LCO for $C<-1/4.$

6. $-(xy^{\prime}(x))^{\prime}-x^{-1}y(x)=\lambda y(x)\;$for all$\;x\in
(0,\infty)$ is LCO at $0$ and LP $+\infty$ (take $\lambda=0$ in $(\ast)$ with
solutions $\cos(\ln(x))$ and $\sin(\ln(x))$. See xamples.f, \#7 (BEZ).

7. $-(xy^{\prime}(x))^{\prime}-xy(x)=\lambda x^{-1}y(x)\;$for all$\;x\in
(0,+\infty)$ is LP at $0$ and LCO at $+\infty$ (take $\lambda=-1/4$ in
$(\ast)$ with solutions $x^{-1/2}\cos(x)$ and $x^{-1/2}\sin(x)$). See
xamples.f, \#6 (Sears-Titchmarsh).

8. $-y^{\prime\prime}(x)+x\sin(x)y(x)=\lambda y(x)$ on $(0,+\infty)$ is R at
$0$ and LP at $+\infty.$ See xamples.f \#30 (Littlewood-McLeod).

\medskip

H5: DEFAULT entry.

The complete range of problems for which SLEIGN2 is applicable can only be
reached by appropriate entries under endpoint classification and boundary
conditions. However, there is a DEFAULT application which requires no detailed
entry of endpoint classification or boundary conditions, subject to:

1) The DEFAULT application CANNOT be used at a LCO endpoint.

2) If an endpoint $a$ is R, then the Dirichlet boundary condition $y(a)=0$ is
automatically used.

3) If an endpoint $a$ is WR, then the following boundary condition is
automatically applied:

(i) if $p(a)=0$, and both $q$ and $w$ are bounded near $a$, then the Dirichlet
boundary condition $y(a)=0$ is used

(ii) if $p(a)>0$, and $q$ and/or $w$ are not bounded near $a$, then the
Neumann boundary condition $y^{\prime}(a)=0$ is used.

If $p(a)=0$, and $q$ and/or $w$ are not bounded near $a$, then no reliable
information, in general, can be given on the DEFAULT boundary condition.

4) If an endpoint is LCNO, then in most cases the principal or Friedrichs
boundary condition is applied (see H6).

5) If an endpoint is LP, then the normal LP procedure is applied (see H7(1.)).

If you choose the DEFAULT condition, then no entry is required for the $u,v$
and $U,V$ boundary condition functions.

\medskip

H6: Limit-circle (LC) boundary conditions.

At an endpoint $a$, the LC type separated boundary condition is of the form
(similar remarks throughout apply to the endpoint $b$ with $U,V$ being
boundary condition functions at $b$)%
\begin{equation}
A1[y,u](a)+A2[y,v](a)=0 \tag{**}%
\end{equation}

\noindent where $y$ is a solution of the differential equation%
\begin{equation}
-(py^{\prime})^{\prime}+qy=\lambda wy\;\text{on}\;(a,b) \tag{*}%
\end{equation}

Here $A1,A2$ are real numbers, not both zero; $u$ and $v$ are boundary
condition functions at $a$; and for real-valued $y$ and $u$ the form
$[y,u](\cdot)$ is defined by%
\[
\lbrack y,u](x)=y(x)(pu^{\prime})(x)-(py^{\prime})(x)u(x)\;\text{for
all}\;x\in(a,b).
\]

If neither endpoint is LP then there are also self-adjoint coupled boundary
conditions. These have a canonical form given by%
\[
Y(b)=\exp(i\alpha)\mathbf{K}Y(a)
\]

\noindent where

(i) $\mathbf{K}$ is a real $2\times2$ matrix with $\det(\mathbf{K})=1$

(ii) the parameter $\alpha$ is restricted to $\alpha\in(-\pi,\pi]$

(iii) $Y$ is the solution column vector $Y(a)=[y(a),(py^{\prime})(a)]^{T}$ at
a regular R endpoint $a$, and $Y$ is the ``singular solution vector''
$Y(a)=[[y,u](a),[y,v](a)]^{T}$ at a singular LC endpoint $a$. Similarly at the
right endpoint $b$ with $U,V.$

The object of this section is to provide help in choosing appropriate
functions $u$ and $v$ in $(\ast\ast)$ (or in choosing $U,V$) given the
differential equation $(\ast)$. Full details of the boundary conditions for
$(\ast)$ are discussed in H7; here it is sufficient to say that the
limit-circle type boundary condition $(\ast\ast)$ must be applied at any
endpoint in the LCNO, LCO classification, but can also be used in the R, WR
classification subject to the appropriate choice of $u$ and $v,$ and $U$ and $V.$

Let $(\ast)$ be R, WR, LCNO, or LCO at endpoint $a$ and choose $c$ in
$(a,b).$. Then \textit{either} $u$ and $v$ are a pair of linearly independent
real solutions of $(\ast)$ on $(a,c]$ for any chosen real value of $\lambda$,
\textit{or }$u$ and $v$ are a pair of real-valued maximal domain functions
defined on $(a,c]$ satisfying $[u,v](a)\neq0.$

The maximal domain $D(a,b)$ is defined by%
\[%
\begin{array}
[c]{lll}%
D(a,b)= & \{f:(a,b)\rightarrow\mathbb{C}: & (i)\;f\text{ and }pf^{\prime}\in
AC_{\text{loc}}(a,b)\\
&  & (ii)\;f\text{ and }w^{-1}(-(pf^{\prime})^{\prime}+qf)\in L^{2}%
((a,b);w)\}.
\end{array}
\]
The domains $D(a,c]$ and $D[c,b)$ are the restrictions of the functions
$D(a,b)$ to the sub-intervals. It is known that for all $f,g\in D(a,c]$ the
limit%
\[
\lbrack f,g](a)=\lim_{x\rightarrow a}[f,g](x)
\]

\noindent exists and is finite. If $(\ast)$ is LCNO or LCO at $a$, then all
solutions of (*) belong to $D(a,c]$ for all values of $\lambda.$ The boundary
condition $(\ast\ast)$ is essential in the LCNO and LCO cases but can also be
used with advantage in some R and WR cases.

In the R, WR, and LCNO cases, but not in the LCO case, the boundary condition
functions can always be chosen so that%
\[
\lim_{x\rightarrow a}\frac{u(x)}{v(x)}=0
\]
and it is recommended that this normalisation be effected, but this is not
essential; this normalization has been entered in the examples given below. In
this case, the boundary condition $[y,u](a)=0$ (\textit{i.e. }$A1=1,A2=0$ in
$(\ast\ast)$) is called the principal or Friedrichs boundary condition at $a$.

In the case when endpoints $a$ and $b$ are, independently, in the R, WR, LCNO,
or LCO classification, it may be that symmetry or other reasons permit one set
of boundary condition functions to be used at both end-points (see xamples.f,
\#1 (Legendre)). In other cases, different pairs must be chosen for each
endpoint (see xamples.f: \#16 (Jacobi), \#18 (Dunsch), and \#19 (Donsch)).

Note that a solution pair $u,v$ is always a maximal domain pair, but not
necessarily vice versa. EXAMPLES:

1. $-y^{\prime\prime}=\lambda y$ on $[0,\pi]$ is R at $0$ and R at $\pi.$ At
$0$, with $\lambda=0$, a solution pair is
\[
u(x)=x,v(x)=1\;\text{for all}\;x\in\lbrack0,\pi].
\]
At $\pi,$with $\lambda=1$, a solution pair is%
\[
U(x)=\sin(x),V(x)=\cos(x)\;\text{for all}\;x\in\lbrack0,\pi].
\]

2. $-(x^{1/2}y^{\prime}(x))^{\prime}=\lambda x^{-1/2}y(x)$ on $(0,1]$ is WR at
$0$ and R at $1.$ (The general solutions of this equation are $u(x)=\cos
(2x^{1/2}\sqrt{\lambda})$ and $v(x)=\sin(2x^{1/2}\sqrt{\lambda})$.)

At $0$, $\lambda=0$, a solution pair is%
\[
u(x)=2x^{1/2},v(x)=1.
\]
At $1$, with $\lambda=\pi^{2}/4$, a solution pair is%
\[
U(x)=\sin(\pi x^{1/2}),V(x)=\cos(\pi x^{1/2}).
\]
Also at 1, with $\lambda=0$, a solution pair is%
\[
U(x)=2(1-x^{1/2}),V(x)=1.
\]

See also xamples.f, \#10 (Weakly Regular).

3. $-((1-x^{2})y^{\prime}(x))^{\prime}=\lambda y(x)$ on $(-1,+1)$ is LCNO at
both endpoints.

At both $\pm1$, $\lambda=0$, a solution pair is%
\[
u(x)=1,v(x)=\ln((1+x)/(1-x))/2
\]

At $+1$, a maximal domain pair is $U(x)=1,V(x)=\ln(1-x)$. At $-1$, a maximal
domain pair is $u(x)=1,v(x)=\ln(1+x).$

See also xamples.f, \#1 (Legendre).

4. $-y^{\prime\prime}(x)-(4x^{2})^{-1}y(x)=\lambda y(x)$ on $(0,+\infty)$ is
LCNO at $0$ and LP at $+\infty$.

At $0$, a maximal domain pair is%
\[
u(x)=x^{1/2},v(x)=x^{1/2}\ln(x).
\]

See also xamples.f, \#2 (Bessel).

5. $-y^{\prime\prime}(x)-5(4x^{2})^{-1}y(x)=\lambda y(x)$ on $(0,+\infty)$ LCO
at $0$ and LP at $+\infty.$

At $0$, $\lambda=0$, a solution pair is%
\[
u(x)=x^{1/2}\cos(\ln(x)),v(x)=\sin(\ln(x)).
\]

See also xamples.f, \#20 (Krall).

6. $-y^{\prime\prime}(x)=\lambda y(x)$ on $(0,+\infty)$ is LCNO at 0 and LP at +infinity.

At $0$, a maximal domain pair is%
\[
u(x)=x,v(x)=1-x\ln(x).
\]

See also xamples.f, \#4(Boyd).

7. $-(x^{-1}y^{\prime}(x))^{\prime}+(kx^{-2}+k^{2}x^{-1})y(x)=\lambda y(x)$ on
$(0,1]$ with $k$ real and $k\neq0$, is LCNO at 0 and R at 1.

At $0$, a maximal domain pair is%
\[
u(x)=x^{2},v(x)=x-k^{-1}%
\]

See also xamples.f, \#8 (Laplace Tidal Wave).

\medskip

H7: General self-adjoint boundary conditions.

Boundary conditions for Sturm-Liouville boundary value problems%
\begin{equation}
-(py^{\prime})^{\prime}+qy=\lambda wy\;\text{on}\;(a,b) \tag{*}%
\end{equation}

\noindent are \textit{either }SEPARATED, with at most one condition at
endpoint $a$ and at most one condition at endpoint $b$,

\noindent\textit{or }COUPLED, when both $a$ and $b$ are, independently, in one
of the end-point classifications R, WR, LCNO, LCO, in which case two
independent boundary conditions are required which link the solution values
near $a$ to those near $b$. The SLEIGN2 program allows for all self-adjoint
boundary conditions; separated self-adjoint conditions and all cases of
coupled self-adjoint conditions.

Separated Conditions: the boundary conditions to be selected depend upon the
classification of the differential equation at the endpoint, say, $a$:

1. If the endpoint a is LP, then no boundary condition is required or allowed.

2. If the endpoint $a$ is R or WR, then a separated self-adjoint boundary
condition is of the form%
\[
A1y(a)+A2(py^{\prime})(a)=0
\]

\noindent where $A1,A2$ are real constants the user must choose, not both zero.

3. If the endpoint a is LCNO or LCO, then a separated boundary condition is of
the form%
\[
A1[y,u](a)+A2[y,v](a)=0
\]

\noindent where $A1,A2$ are real constants the user must choose, not both
zero; here $u,v$ are the pair of boundary condition functions the user has
previously selected when the input Fortran file was being prepared with makepqw.f.

4. If the endpoint a is LCNO and the boundary condition pair $u,v$ has been
chosen so that%
\[
\lim_{x\rightarrow a}\frac{u(x)}{v(x)}=0
\]

\noindent(which is always possible), then $A1=1,A2=0$ (\textit{i.e.
}$[y,u](a)=0$) gives the principal (Friedrichs) boundary condition at $a.$

5. If $a$ is R or WR and boundary condition functions $u,v$ have been entered
in the Fortran input file, then 3 and 4 above apply to entering separated
boundary conditions at such an endpoint; the boundary conditions in this form
are equivalent to the point-wise conditions in 2 (subject to care in choosing
$A1,A2$). This singular form of a regular boundary condition may be
particularly effective in the WR case if the boundary condition form in 2
leads to numerical difficulties. Conditions 2,3,4,5 apply similarly at
endpoint $b$ (with $U,V$ as the boundary condition functions at $b$).

6. If $a$ is R, WR, LCNO, or LCO and $b$ is LP, then only a separated
condition at $a$ is required and allowed (or instead at $b$ if $a$ and $b$ are interchanged).

7. If both endpoints $a$ and $b$ are LP, then no boundary conditions are
required or allowed.

The indexing of eigenvalues for boundary value problems with separated
conditions is discussed in H13.

Coupled Conditions:

8. Coupled regular self-adjoint boundary conditions on $(a,b)$ apply only when
both endpoints $a$ and $b$ are R or WR.

\medskip

H8: Recording the results.

If you choose to have a record kept of the results, then the following
information is stored in a file with the name you select:

1. The file name.

2. The header line prompted for (up to 32 characters of your choice).

3. The interval $(a,b)$ selected by the user.

For SEPARATED boundary conditions:

4. The endpoint classification.

5. A summary of coefficient information at WR, LCNO, LCO endpoints.

6. The boundary condition constants $(A1,A2),(B1,B2)$ if entered.

7. (NUMEIG,EIG,TOL) or (NUMEIG1,NUMEIG2,TOL), as entered.

For COUPLED boundary conditions:

8. The boundary condition parameter $\alpha$ and the coupling matrix
$\mathbf{K}$, see H6.

For ALL self-adjoint boundary conditions:

9. The computed eigenvalue, EIG, and its estimated accuracy, TOL.

10. IFLAG reported (see H15).

\medskip

H9: Type and choice of interval.

You may enter any interval $(a,b)$ for which the coefficients $p,q,w$ are well
defined by your Fortran statements in the input file, provided that $(a,b)$
contains no interior singularities.

\medskip

H10: Entry of endpoints.

Endpoints a and b should generally be entered as real numbers to an
appropriate number of decimal places.

\medskip

H11: endpoint values of $p,q,w.$

The program SLEIGN2 needs to know whether the coefficient functions $p,q,w$ as
defined by the Fortran expressions entered in the input file, can be evaluated
numerically without running into difficulty. If, for example, either $q$ or
$w$ is unbounded at $a$, or $p(a)=0$, then SLEIGN2 needs to know this
information so that $a$ is not chosen for functional evaluation.

\medskip

H12: Initial value problems.

The initial value problem facility for Sturm-Liouville problems%
\begin{equation}
-(py^{\prime})^{\prime}+qy=\lambda wy\;\text{on}\;(a,b) \tag{*}%
\end{equation}

\noindent allows for the computation of a solution of $(\ast)$ with a
user-chosen value $\lambda$ and any one of the following initial conditions:

1. From endpoint $a$ of any classification except LP towards endpoint $b$ of
any classification

2. From endpoint $b$ of any classification except LP back towards endpoint $a$
of any classification

3. From endpoints $a$ and $b$ of any classifications except LP towards an
interior point of $(a,b)$ selected by the program.

Initial values at $a$ are of the form $y(a)=\alpha1,(py^{\prime})(a)=\alpha2$
when $a$ is R or WR; and $[y,u](a)=\alpha1,[y,v](a)=\alpha2$ when $a$ is LCNO
or LCO.

Initial values at $b$ are of the form $y(b)=\beta1,(py^{\prime})(b)=\beta2$
when $b$ is R or WR; and $[y,u](b)=\beta1,[y,v](b)=\beta2$ when $b$ is LCNO or LCO.

In $(\ast)$, $\lambda$ is a user-chosen real number; while in the above
initial values, $(\alpha1,\alpha2)$ and $(\beta1,\beta2)$ are user-chosen
pairs of real numbers not both zero.

In the initial value case 3 above when the interval $(a,b)$ is finite, the
interior point selected by the program is generally near the midpoint of
$(a,b)$; when $(a,b)$ is infinite, no general rule can be given. Also if,
given $(\alpha1,\alpha2)$ and $(\beta1,\beta2)$, the $\lambda$ chosen is an
eigenvalue of the associated boundary value problem, the computed solution may
not be the corresponding eigenfunction -- the signs of the computed solutions
on either side of the interior point may be opposite.

The output for a solution of an initial value problem is in the form of stored
numerical data which can be plotted on the screen (see H16), or printed out in
graphical form if graphics software is available.

\medskip

H13: Indexing of eigenvalues.

The indexing of eigenvalues is an automatic facility in SLEIGN2. The following
general results hold for the separated boundary condition problem (see H7):

1. If neither endpoint $a$ or $b$ is LP or LCO, then the spectrum of the
eigenvalue problem is discrete (eigenvalues only), simple (eigenvalues all of
multiplicity 1), and bounded below with a single cluster point at $+\infty.$
The eigenvalues are indexed as $\{\lambda_{n}:n=0,1,2,...\}$, where
$\lambda_{n}<\lambda_{n+1}$ $(n=0,1,2,..)$, $\lim_{n\rightarrow+\infty}%
\lambda_{n}=+\infty$; and if $\{\psi_{n}:n=0,1,2,...\}$ are the corresponding
eigenfunctions, then $\psi_{n}$ has exactly $n$ zeros in the open interval $(a,b).$

2. If neither endpoint $a$ or $b$ is LP but at least one endpoint is LCO, then
the spectrum is discrete and simple as for 1, but with cluster points at both
$\pm\infty.$ The eigenvalues are indexed as $\{\lambda_{n}:n=0,\pm
1,\pm2,...\}$, where $\lambda_{n}<\lambda_{n+1}$ $($for $n<n+1$ and
$n=0,\pm1,\pm2,..),$ with $\lambda_{0}$ the smallest non-negative eigenvalue;
$\lim_{n\rightarrow-\infty}\lambda_{n}=-\infty$ and $\lim_{n\rightarrow
+\infty}\lambda_{n}=+\infty$; and if $\{\psi_{n}:n=0,\pm1,\pm2,...\}$ are the
corresponding eigenfunctions, then every $\psi_{n}$ has infinitely many zeros
in $(a,b).$

3. If one or both endpoints is LP, then there can be one or more intervals of
continuous spectrum for the boundary value problem in addition to some
(necessarily simple) eigenvalues. For these essentially more difficult
problems, SLEIGN2 can be used as an investigative tool to give qualitative and
possibly quantitative information on the spectrum.

For example, if a problem has continuous spectrum starting at a real number
$L$, then there may be no eigenvalues below $L$, any finite number of
eigenvalues below $L$, or an infinite (but countable) number of eigenvalues
below $L.$ SLEIGN2 can be used to compute $L$ (see the paper BEWZ on the
SLEIGN2 home page for an algorithm to compute $L$), and to determine the
number of these eigenvalues and compute them. In this respect, see xamples.f:
\#13 (Hydrogen Atom), \#17 (Morse Oscillator), \#21 (Fourier), and \#27
(J\"{o}rgens) as examples of success; and \#2 (Mathieu), \#14 (Marletta), and
\#28 (Behnke-Goerisch) as examples of failure.

The problem need not have a continuous spectrum, in which case if its discrete
spectrum is bounded below, then the eigenvalues are indexed and the
eigenfunctions have zero counts as in 1. If, on the other hand, the discrete
spectrum is unbounded below, then all the eigenfunctions have infinitely many
zeros in the open interval $(a,b)$. SLEIGN2 can, in principle, compute these
eigenvalues if neither endpoint is LP although this is a computationally
difficult problem. Note however, that SLEIGN2 has no algorithm when the
spectrum is discrete, unbounded above and below and one endpoint is LP, as in
xamples.f \#30.

In respect to the five classes of endpoints, the following identified examples
from xamples.f illustrate the spectral property of these boundary value problems:

1. Neither endpoint is LP or LCO:

\#1 (Legendre), \#2 (Bessel) with $-1/4<C<3/4,$ \#4 (Boyd), \#5 (Latzko).

2. Neither endpoint is LP, but at least one is LCO:

\#6 (Sears-Titchmarsh), \#7 (BEZ), \#19 (Donsch)

3. At least one endpoint is LP:

\#13 (Hydrogen Atom), \#14 (Marletta), \#20 (Krall), \#21 (Fourier) on [0,infinity).

\medskip

H14: Entry of eigenvalue index, initial guess, and tolerance.

For all self-adjoint boundary condition problems (see H7), SLEIGN2 calls for
input information options to compute \textit{either}

1. a single eigenvalue, \textit{or}

2. a series of eigenvalues.

\noindent In each case indexing of eigenvalues is called for (see H13).

1 above asks for data triples NUMEIG, EIG, TOL separated by commas. Here
NUMEIG is the integer index of the desired eigenvalue; NUMEIG can be negative
only when the problem is LCO at one or both endpoints. EIG allows for the
entry of an initial guess for the requested eigenvalue (if an especially good
one is available), or can be set to 0 in which case an initial guess is
generated by SLEIGN2 itself. TOL is the desired accuracy of the computed
eigenvalue. It is an absolute accuracy if the magnitude of the eigenvalue is 1
or less, and is a relative accuracy otherwise. Typical values for TOL might be
0.001 for moderate accuracy and 0.0000001 for high accuracy in single
precision. If TOL is set to 0, the maximum achievable accuracy is requested.

If the input data list is truncated with a ``\thinspace/\thinspace'' after
NUMEIG or EIG, then the remaining elements default to 0.

2 above asks for data triples NUMEIG1, NUMEIG2, TOL separated by commas. Here
NUMEIG1 and NUMEIG2 are the first and last integer indices of the sequence of
desired eigenvalues, with NUMEIG1
%TCIMACRO{\TEXTsymbol{<}}%
%BeginExpansion
$<$%
%EndExpansion
NUMEIG2; they can be negative only when the problem is LCO at one or both
endpoints. TOL is the desired accuracy of the computed eigenvalues. It is an
absolute accuracy if the magnitude of an eigenvalue is 1 or less, and is a
relative accuracy otherwise. Typical values for TOL might be 0.001 for
moderate accuracy and 0.0000001 for high accuracy in single precision. If TOL
is set to 0, the maximum achievable accuracy is requested.

If the input data list is truncated with a ``\thinspace/\thinspace'' after
NUMEIG2, then TOL defaults to 0.

For COUPLED self-adjoint boundary condition problems (see H7 and H17), SLEIGN2
also reports which eigenvalues are double. Double eigenvalues can occur only
for coupled boundary conditions with the parameter $\alpha=0$ or $\pi.$

\medskip

H15: IFLAG information.

All results are reported by SLEIGN2 with a flag identification. There are four
values of IFLAG:

1. The computed eigenvalue has an estimated accuracy within the

tolerance requested.

2. The computed eigenvalue does not have an estimated accuracy within the
tolerance requested, but is the best the program could obtain.

3. There seems to be no eigenvalue of index equal to NUMEIG.

4. The program has been unable to compute the requested eigenvalue.

\medskip

H16: Plotting.

After computing a single eigenvalue (see H14(1.)), but not a sequence of
eigenvalues (see H14(2.)), the eigenfunction can be plotted for separated
conditions and for coupled ones with $\alpha=0,\pi.$ If this data is desired,
respond $y$ when asked and SLEIGN2 will then compute eigenfunction data and
store them for subsequent use.

The user can ask that the eigenfunction data be in the form of either points
$(x,y)$ for $x\in(a,b)$, or points $(t,y)$ for $t$ in the standardized
interval $(-1,+1)$ mapped onto from $(a,b)$; the $t$-choice can be especially
helpful when the original interval is infinite. Additionally, the user can ask
for a plot of the so-called Pr\"{u}fer angle, in the $x$- or $t$-variables.

In both forms, once the choice has been made of the function to be plotted, a
crude plot is displayed on the monitor screen and the user is asked whether or
not to save the computed plot points in a file.

\medskip

H17: Indexing of eigenvalues for coupled self-adjoint problems.

The indexing of eigenvalues is an automatic facility in SLEIGN2. The following
general result holds for coupled boundary condition problems (see H7):

The spectrum of the eigenvalue problem is discrete (eigenvalues only). In
general the spectrum is not simple, but no eigenvalue exceeds multiplicity
$2$. The eigenvalues are indexed as $\{\lambda_{n}:n=0,1,2,...\}$ where
$\lambda_{n}<\lambda_{n+1}$ for $n=0,1,2,...,$and $\lim_{n\rightarrow+\infty
}\lambda_{n}=+\infty$ if neither endpoint is LCO. If one or both endpoints are
LCO the eigenvalues cluster at both $-\infty$ and at $+\infty,$ and all
eigenfunctions have infinitely many zeros.

If neither endpoint is LCO and $\alpha=0,\pi$, then the $n$-th eigenfunction
has $n-1,n,n+1$ zeros in the half open interval $[a,b)$ (also in $(a,b]$). All
three possibilities occur. Recall that in the case of double eigenvalues,
although the $n$-th eigenvalue is well defined there is an ambiguity about
which solution is declared as the $n$-th eigenfunction. If $\alpha\neq0,\pi$,
then the eigenfunction is non-real and has no zero in $(a,b)$, but each of the
real and imaginary parts of the $n$-th eigenfunction have the same zero
properties as mentioned above when $\alpha=0,\pi.$

The following identified examples from xamples.f are of special interest:

\#11 (Plum) on $[0,\pi],$ \#21 (Fourier) on $[0,\pi],$ \#25 (Meissner) on $[-1/2,+1/2].$
\end{document}
SHAR_EOF
fi # end of overwriting check
if test -f 'intro.tex'
then
	echo shar: will not over-write existing file "'intro.tex'"
else
cat << "SHAR_EOF" > 'intro.tex'
%% This document created by Scientific Word (R) Version 3.5

\documentclass[12pt]{amsart}%
\usepackage{graphicx}
\usepackage{amscd}
\usepackage{amsmath}%
\usepackage{amsfonts}%
\usepackage{amssymb}
%TCIDATA{OutputFilter=latex2.dll}
%TCIDATA{CSTFile=amsartci.cst}
%TCIDATA{Created=Mon Sep 06 21:20:17 1999}
%TCIDATA{LastRevised=Sunday, April 01, 2001 11:12:39}
%TCIDATA{<META NAME="GraphicsSave" CONTENT="32">}
%TCIDATA{<META NAME="DocumentShell" CONTENT="Journal Articles\amsnumb">}
%TCIDATA{Language=British English}
\newtheorem{theorem}{Theorem}
\theoremstyle{plain}
\newtheorem{acknowledgement}{Acknowledgement}
\newtheorem{algorithm}{Algorithm}[section]
\newtheorem{axiom}{Axiom}[section]
\newtheorem{case}{Case}[section]
\newtheorem{claim}{Claim}[section]
\newtheorem{conclusion}{Conclusion}[section]
\newtheorem{condition}{Condition}[section]
\newtheorem{conjecture}{Conjecture}[section]
\newtheorem{corollary}{Corollary}[section]
\newtheorem{criterion}{Criterion}[section]
\newtheorem{definition}{Definition}[section]
\newtheorem{example}{Example}[section]
\newtheorem{exercise}{Exercise}[section]
\newtheorem{lemma}{Lemma}[section]
\newtheorem{notation}{Notation}[section]
\newtheorem{problem}{Problem}[section]
\newtheorem{proposition}{Proposition}[section]
\newtheorem{remark}{Remark}[section]
\newtheorem{solution}{Solution}[section]
\newtheorem{summary}{Summary}[section]
\numberwithin{equation}{section}
\numberwithin{theorem}{section}
\newcommand{\thmref}[1]{Theorem~\ref{#1}}
\newcommand{\secref}[1]{\S\ref{#1}}
\newcommand{\lemref}[1]{Lemma~\ref{#1}}
\setlength{\oddsidemargin}{0.0in}
\setlength{\evensidemargin}{0.0in}
\setlength{\textwidth}{6.5in}
\setlength{\textheight}{9.0in}
\setlength{\headsep}{0.25in}
\setlength{\headheight}{0.0in}


\begin{document}
\title[Introduction to SLEIGN2]{Introduction to SLEIGN2}
\author{P.B. Bailey}
\author{W.N. Everitt}
\author{A. Zettl}
\address{Department of Mathematical Sciences, Northern Illinois University, DeKalb, IL
60115-2888, USA}
\date{01 March 2001 (File: intro.tex).}
\maketitle


The main purpose of this code is to compute eigenvalues and eigenfunctions of
regular and singular self-adjoint Sturm-Liouville problems (SLP) and to
approximate the continuous spectrum in the singular case. For a general
description of the analytical and numerical properties of the SLEIGN2 code see
\cite{BEZ} and \cite{BEZ3}.

These problems consist of a second order linear differential equation
\[
-(py^{\prime})^{\prime}+qy=\lambda wy\;\text{on}\;(a,b)
\]
together with boundary conditions (BC). The nature of the BC depends on the
regular or singular classification of the end points a and b. For both cases
the BC fall into two major classes: separated and coupled. The former are two
separate conditions, one at each endpoint; the latter are two coupled
conditions linking the values of the solution near the two endpoints, e.g.
periodic and semi-periodic boundary conditions. SLEIGN2 seems to be the only
general purpose code available for arbitrary self-adjoint BC, separated or
coupled, and for both regular and singular problems.

A number $\lambda$ for which there is a nontrivial solution satisfying the BC
is called an eigenvalue and such a solution is a (corresponding)
eigenfunction. If one or both endpoints are LP (see below or section 4 of HELP
for a definition) there may be points $\lambda$ in the spectrum in addition to
eigenvalues i.e. there may be continuous spectrum.

In the theory of SLP the coefficients $1/p$ and $q$ and the weight function
$w$ are assumed to be real-valued and locally Lebesgue integrable on $(a,b).$

To meet the needs of numerical computing techniques we make the stronger assumptions:

(i) The interval $(a,b)$ of R may be bounded or unbounded

(ii) $p,q$ and $w$ are real-valued functions on $(a,b)$

(iii) $p,q$ and $w$ are piecewise continuous on $(a,b)$

(iv) $p$ and $w$ are strictly positive on $(a,b)$.

For better error analysis in the numerical procedures, condition (iii) above
is replaced with

(iii)$^{\prime}$ $p,q$ and $w$ are four times continuously differentiable on
$(a,b)$.

To study a SLP using operator theory one associates a self-adjoint operator in
the weighted Hilbert space of square-integrable functions, with respect to the
weight $w$, on $(a,b)$, with each SLP in such a way that the spectrum of the
problem is the spectrum of the operator. In the case of a regular problem the
spectrum consists entirely of eigenvalues and these are bounded below (when
$p>0$ and $w>0$). This is still so for the case when each endpoint is either
regular (R) or singular limit-circle nonoscillatory (LCNO). In case one
endpoint is limit-circle oscillatory (LCO) and the other is not limit-point
(LP) then there are still only eigenvalues in the spectrum but these are not
bounded below. (The spectrum is never bounded above.)

If one or both endpoints is LP the spectrum may be extremely complicated.
There may be no eigenvalues, finitely many, or infinitely many. Some may be
embedded in the continuous spectrum. For $p=1,w=1,q(x)=\sin(x)$ on
$(-\infty,+\infty)$ there are no eigenvalues and the continuous spectrum
consists of the union of an infinite number of disjoint compact intervals.
(SLEIGN2 can be used to approximate this spectrum - see example 12 in the file
xamples.tex and the references quoted there.)

See the HELP file help.tex for a definition of the terms regular (R),
limit-circle (LC), limit-circle non-oscillatory (LCNO), limit-circle
oscillatory (LCO), limit-point (LP).

SLP problems are classified into various classes based on the classification
of the endpoints and on whether the boundary conditions are separated (S) or
coupled (C). We have the following categories:

1. R/R, S

2. R/R, C

3. R/LCNO or LCNO/R, S

4. R/LCNO or LCNO/R, C

5. R/LCO or LCO/R , S

6. R/LCO or LCO/R , C

7. LCNO/LCO or LCO/LCNO or LCO/LCO, S

8. LCNO/LCO or LCO/LCNO or LCO/LCO, C

9. LP/R or LP/LCNO or LP/LCO or R/LP or LCNO/LP or LCO/LP

10. LP/LP

For 9 there is only a separated condition at the non-LP endpoint and for 10
there are no boundary conditions at either end.

There are only two other major general purpose codes for computing eigenvalues
and eigenfunctions of Sturm-Liouville problems : SLEDGE (Fulton and Pruess)
and the earlier code SLEIGN (Bailey \textit{et al).} (There was another code,
written by Pryce and Marletta, available from the NAG library but the current
NAG library no longer includes such a code.)

SLEDGE uses a method based on piecewise constant approximations of the
coefficients of the differential equation; SLEIGN and SLEIGN2 both use the
Pr\"{u}fer transformation.

Both SLEIGN and SLEDGE are designed for separated regular boundary conditions
and both have a mechanism to automatically handle endpoints which are either
regular or singular but non-oscillatory. In the latter case if an endpoint is
LCNO the Friedrichs condition is usually, but not always, the one chosen by
the code. SLEDGE can also determine the LP/LC classification but only in a
restricted number of cases (the method requires that the coefficients $p,q$
and $w$ are analytic on $(a,b)$ and depends on the Frobenius method for series
solutions at regular singular endpoints); SLEIGN and SLEIGN2 do not have such
a facility but, for the sake of generality, rely on the user input of this information.

SLEIGN2 is the only general purpose code in existence which can, in principle,
handle arbitrary self-adjoint, separated or coupled, regular or singular,
boundary conditions, see the HELP file help.tex and \cite{BEZ3}. Problems with
coupled boundary conditions, see\cite{BEZ2}, at singular endpoints are
difficult to handle numerically, especially if one or both of the endpoints is LCO.

In addition to the above mentioned capabilities to compute eigenvalues and
eigenfunctions SLEIGN2 also computes the solution of an initial value problem
with the users choice of $\lambda$ and either a regular or singular initial condition.

When combined with the algorithm established in \cite{BEWZ}, SLEIGN2 can be
used to approximate the continuous spectrum.

An important feature of the SLEIGN2 program is its user friendly interface
contained in makepqw.f and drive.f. Users who want to bypass this interface
and use their own driver may wish to look at a sample driver; two such drivers
are provided, see (7) and (8) below.

The whole package consists of the following files:

1. A brief ``readme'' file readme.txt with basic information on how to run the code.

2. This introduction file intro.tex.

3. makepqw.f - This is an interactive Fortran file to input the coefficient
functions $p,q,w$ and, if necessary, the functions $u,v$ which define the
singular boundary conditions

4. drive.f - This is an interactive Fortran file containing the driver, a HELP
file help.tex, and a ``user friendly'' interface.

5. sleign2.f - The main code for the computation of eigenvalues and eigenfunctions.

6. xamples.f - A Fortran file with 32 examples ready to run. These examples
were chosen to illustrate various features of the code.

7. sepdr.f - A sample driver for separated regular or singular boundary conditions.

8. coupdr.f - A sample driver for coupled regular or singular boundary conditions.

9. xamples.tex - A LaTeX file containing information about the 32 examples.

10. help.tex - A LaTeX file with information about endpoint classifications,
boundary conditions etc. It is a separate text file and it can also be
accessed from both makepqw.f and drive.f.

11. autoinput.txt- A file describing an ``automatic'' method for using SLEIGN2
which avoids the user friendly ``question and answer'' format; this is
recommended for experienced users only.

When an eigenfunction has been computed it is stored. If it is real-valued, it
can be examined:

(i) by printing out the numerical data

(ii) by using the discrete graph plotter in the program

(iii) by using a local graph plotter.

Full details are given in HELP; see the file help.tex.

All six of the Fortran files in the SLEIGN2 package are in single precision.

\begin{thebibliography}{9}                                                                                                %

\bibitem {BEZ}P.B. Bailey, W.N. Everitt and A. Zettl, \emph{Computing
eigenvalues of singular Sturm-Liouville problems}, Results in Mathematics,
\textbf{20 }(1991), 391-423.

\bibitem {BEZ2}P.B. Bailey, W.N. Everitt and A. Zettl, \emph{Regular and
singular Sturm-Liouville problems with coupled boundary conditions}, Proc.
Royal Soc. Edinburgh (A) \textbf{126} (1996), 505-514.

\bibitem {BEZ3}P.B. Bailey, W.N. Everitt and A. Zettl, \emph{The SLEIGN2
Sturm-Liouville code.} (To appear in ACM Trans. Math. Software).

\bibitem {BEWZ}P.B. Bailey, W.N. Everitt, J. Weidmann and A. Zettl,
\emph{Regular approximation of singular Sturm-Liouville problems}. Results in
Mathematics \textbf{23} (1993), 3-22.
\end{thebibliography}
\end{document}
SHAR_EOF
fi # end of overwriting check
if test -f 'readme.txt'
then
	echo shar: will not over-write existing file "'readme.txt'"
else
cat << "SHAR_EOF" > 'readme.txt'

 
                   SLEIGN2: THE README.TXT FILE
  
   01 March 2001; P.B. BAILEY, W.N. EVERITT AND A. ZETTL
  
  PLEASE read carefully this readme.txt file and the intro.tex file before 
 using the SLEIGN2 package for the first time.
  
 	There are eleven files in the SLEIGN2 package as follows:
 
 	Two ASCII files:
 	autoinput.txt  readme.txt
 
 	Six FORTRAN files:
 	coupdr.f  drive.f  makepqw.f  sepdr.f  sleign2.f  xamples.f
 
 	Three AMS-LaTeX files (these files can be compiled in the UNIX latex
 compiler, and then printed out in hard copy):
 	help.tex  intro.tex  xamples.tex   
  
   To run one of the examples in the xamples.f file in an UNIX environment 
 with a Fortran77 (or Fortran77) compiler, enter the following command:
  
 f77 xamples.f  drive.f  sleign2.f  -o  xamples.x
 (Replace f77 by f90 if you want to use the Fortran90 compiler.)
  
   This will create the executable file xamples.x and object files drive.o
 and sleign2.o. Now run xamples.x whenever you want to work an example from 
 the list of examples in the xamples.f file. The hard printed copy of the
 file xamples.tex provides detailed information on each one of the 32
 chosen examples in the file xamples.x.
  
   To run your own problem proceed as follows:
  Step 1: Enter the command 
        f77  makepqw.f  -o  makepqw.x
  Step 2: Run
        makepqw.x  
   (This interactive program will ask you for a file name - this must end 
 in .f - for example, if your chosen problem has the name bloggs then enter 
 bloggs.f). This file, after makepqw.x has been run, will contain the 
 subroutines for p,q,w and, if necessary, the functions u, v and U, V which 
 are used to define singular limit-circle boundary conditions at one or both 
 endpoints.
  Step 3: Enter the command
        f77  bloggs.f  drive.f  sleign2.f  -o  bloggs.x
   (drive.f and sleign2.f can be replaced with drive.o and sleign2.o 
 if these .o files are available to speed up the compilation; these -o 
 files are created by the first compilation.)
  Step 4: Run
        bloggs.x  
   (The user is asked to provide information for the code to run
 bloggs.x, for example: boundary conditions, eigenvalue indexes, 
 numerical tolerances, name for report file if desired. See the autoinput.txt 
 file for instructions on how to automate the input of the required 
 information).
  
  The file sepdr.f is a sample driver for separated regular and singular
 boundary conditions; coupdr.f is a sample driver program for coupled regular
 and singular boundary conditions. Experienced users who want to bypass the 
 extensive user friendly interface provided in drive.f and makepqw.f, and
 use their own driver may wish to look at these two sample drivers.

  Note that the above procedures may have to be modified for non UNIX 
 environments, e.g. DOS or APPLE.
 
  Note also that in running xamples.x or bloggs.x the user has access to
 the interactive help device; at any point where the program halts  
 type h <ENTER> for access; to return to the point in the program where
 help was accessed, type r <ENTER>. Additional information on the help 
 device can be found in the intro.tex file. The whole of the help data
 can be printed out separately from the file help.tex and the user is
 advised to have a copy available for consultation when working with the
 code files for the first time.
 
  See the autoinput.txt file for instructions on how to bypass these 
 program halts in xamples.x and bloggs.x
  
  All six of the FORTRAN .f files are supplied in single precision; 
 to convert these files to double precision replace the string `  REAL' by 
 the string `  DOUBLE PRECISION' throughout. (Note the two spaces in front of 
 REAL and in front of DOUBLE PRECISION - these spaces are important.)  
 In UNIX this can be effected within the vi editor as follows:
  
 :1,$ s/  REAL/  DOUBLE PRECISION
  
  It is recommended that the user try the program in single precision and 
 switch to double precision as required.
  
  Additional information on the SLEIGN2 package can be found in the intro.tex
 and help.tex files. See also a hard printed copy of the file bailey.tex,
 available in the Zettl/Sleign2 directory detailed below, which contains
 an account of the analytical and numerical properties of the SLEIGN2
 code, together with an extensive list of references.
  
  All of the eleven files in the SLEIGN2 package can be obtained by anonymous
 ftp from the computer
        ftp.math.niu.edu
 and the directory
        /pub/papers/Zettl/Sleign2
  
  A sample session using the traditional UNIX or DOS ftp client is as 
 follows:
  
  Step 1. At the UNIX or DOS prompt type
              fop  ftp.math.niu.edu
  
  Step 2. In response to the request for user id type
              fop
  
  Step 3. In response to the request for username type your e-mail address;
             for example:
             w.n.everitt@bham.ac.uk
  
  Step 4. At the ftp> prompt enter:
             cd  pub/papers/Zettl/Sleign2
             get  readme.txt
             quit
  
  This will transfer the readme.txt file to your current directory and close 
 the connection.
  
  At the ftp> prompt you can also enter 'ls' to see the listing of other 
 files available in the Sleign2 directory.
  
  Users who prefer World Wide Web software such as Netscape 
 or lynx can specify the URL
            ftp://ftp.math.niu.edu/pub/papers/Zettl/Sleign2
 to access the Sleign2 directory; this software will then show the list 
 of files available in it.
 
  Replacing the directory "Sleign2" with the directory "Pub papers" and
 repeating the above procedures will make available a number of recent papers
 which are related to the SLEIGN2 package. For example, the paper "BEWZ" 
 (Bailey, Everitt, Weidmann and Zettl) contains some results which can be 
 combined with SLEIGN2 to approximate the continuous (essiential) spectrum 
 of singular limit-point Sturm-Liouville problems.
  
  All eleven files of the SLEIGN2 package and a number of recent publications
 related to it, can also be accessed through the web page:
           http://www.math.niu.edu/~zettl/SL2/
  
  All suggestions, comments and criticisms are welcome; please send all 
 comments to Tony Zettl at
          zettl@math.niu.edu
  
     Paul Bailey, Norrie Everitt and Tony Zettl 
 (with the assistance of Burt Garbow.)
  
  Acknowledgement. The authors are grateful to their colleagues Eric Behr, 
 Qingkai Kong and Hongyou Wu for help and advice at a number of stages 
 in the development of this program. Some of the theoretical underpinnings 
 of the algorithm for coupled boundary conditions were obtained jointly 
 with Michael Eastham, Qingkai Kong and Hongyou Wu.

  A special thanks to Eric Behr for his help throughout the development of
 the code, for setting up the public access through the Internet and the
 world wide web, and for informed advice.

 Department of Mathematical Sciences, Northern Illinois University 
 DeKalb, IL 60115-2888, USA



SHAR_EOF
fi # end of overwriting check
if test -f 'xamples.tex'
then
	echo shar: will not over-write existing file "'xamples.tex'"
else
cat << "SHAR_EOF" > 'xamples.tex'
%% This document created by Scientific Word (R) Version 3.5

\documentclass[12pt]{amsart}%
\usepackage{graphicx}
\usepackage{amscd}
\usepackage{amsmath}%
\usepackage{amsfonts}%
\usepackage{amssymb}
%TCIDATA{OutputFilter=latex2.dll}
%TCIDATA{CSTFile=amsartci.cst}
%TCIDATA{Created=Tue Mar 09 09:05:39 1999}
%TCIDATA{LastRevised=Wednesday, April 04, 2001 12:27:55}
%TCIDATA{<META NAME="GraphicsSave" CONTENT="32">}
%TCIDATA{<META NAME="DocumentShell" CONTENT="Journal Articles\amslatex">}
%TCIDATA{Language=British English}
\newtheorem{theorem}{Theorem}
\theoremstyle{plain}
\newtheorem{acknowledgement}{Acknowledgement}
\newtheorem{algorithm}{Algorithm}
\newtheorem{axiom}{Axiom}
\newtheorem{case}{Case}
\newtheorem{claim}{Claim}
\newtheorem{conclusion}{Conclusion}
\newtheorem{condition}{Condition}
\newtheorem{conjecture}{Conjecture}
\newtheorem{corollary}{Corollary}
\newtheorem{criterion}{Criterion}
\newtheorem{definition}{Definition}
\newtheorem{example}{Example}
\newtheorem{exercise}{Exercise}
\newtheorem{lemma}{Lemma}
\newtheorem{notation}{Notation}
\newtheorem{problem}{Problem}
\newtheorem{proposition}{Proposition}
\newtheorem{remark}{Remark}
\newtheorem{solution}{Solution}
\newtheorem{summary}{Summary}
\numberwithin{equation}{section}
\newcommand{\thmref}[1]{Theorem~\ref{#1}}
\newcommand{\secref}[1]{\S\ref{#1}}
\newcommand{\lemref}[1]{Lemma~\ref{#1}}
\setlength{\oddsidemargin}{0.0in}
\setlength{\evensidemargin}{0.0in}
\setlength{\textwidth}{6.5in}
\setlength{\textheight}{8.5in}
\setlength{\headsep}{0.25in}
\setlength{\headheight}{0.0in}


\begin{document}
\title[SLEIGN2]{SLEIGN2\\Commentary on the individual examples in xamples.f}
\author{P.B. Bailey}
\address{P.B. Bailey, c/o Department of Mathematical Sciences, Northern Illinois
University, DeKalb, IL 60155-2888, USA}
\email{70621.3674@compuserve.com}
\author{W.N. Everitt}
\address{W.N. Everitt, School of Mathematics and Statistics, University of Birmingham,
Edgbaston, Birmingham B15 2TT, England, UK}
\email{w.n.everitt@bham.ac.uk}
\author{A. Zettl}
\address{A. Zettl, Department of Mathematical Sciences, Northern Illinois University,
DeKalb, IL 60155-2888, USA}
\email{zettl@math.niu.edu}
\date{01 March 2001 (File: xamples.tex)\quad\AmS  -\LaTeX  {}; prepared in
\textsl{Scientific Word\/}{}}
\maketitle


\section{Introduction}

The examples in this commentary have been chosen to illustrate the
capabilities and limitations of the program SLEIGN2. Many of the examples have
been chosen from special cases of the well known and well studied ``special
functions'' of mathematical analysis. All possible cases of endpoint
classification are represented; all types of self-adjoint boundary conditions
are included, \textit{i.e.} regular or singular, and separated or coupled. In
the limit-circle case examples are given for which the endpoints may be
oscillatory or non-oscillatory.

For a general account of both the analytical and numerical properties of the
SLEIGN2 code see the paper by Bailey, Everitt and Zettl \cite{BEZ3}.

For all 32 examples in this commentary the following data have been entered:

(i) the Sturm-Liouville differential equation and associated interval on the
real line $\mathbb{R}$

(ii) the range of any parameters in the differential equation; this serves to
remind the reader that numerical values for these parameters have to be
entered in some of the examples given in xamples.x, for any such example to run

(iii) the endpoint classification of the differential equation, in the
relevant Hilbert function space $L^{2}((a,b);w)$

(iv) the boundary condition functions $u,v$ required for any LCNO or LCO endpoint

(v) comments on any particular features of the numbered example.

The data in items (i) to (iv) above can also be found in the file xamples.f,
but this search requires scrolling through the file as the data items, for any
particular example, are located in different sections.

For some of these examples it is possible to give explicit information on the
spectrum of associated boundary value problems; this can take the form of
providing explicit formulas for eigenvalues against which the program
calculated results can be compared.

In all cases of limit-circle endpoints, boundary condition functions $u$ and
$v$ have been entered as part of the example data. In the case of limit-circle
non-oscillatory endpoints we use the convention that the boundary condition
function $u$ determines the principal or Friedrichs boundary condition.

On selecting a numbered example in the file xamples.x, the differential
equation is displayed in Fortran, and details of the endpoint classification
given. If information on the form of the boundary condition functions $u$ and
$v$ is required then the user should scroll separately through the file
xamples.f to the appropriate numbered part of the $u,v$ section or refer to
the information in this commentary.

Some regular and weakly regular problems can be more successfully run using
the limit-circle non-oscillatory (LCNO) algorithm; details are given below for
some of the examples.

It should be noted that for limit-circle oscillatory problems it is sometimes
difficult to compute numerically more than a few of the eigenvalues. This is
due, at least in part, to the rapid growth of the eigenvalues in both the
positive and the negative directions; but particularly in the negative direction.

The Laguerre problem, Example 22, has a discrete spectrum and for one
particular boundary condition the eigenvalues are known explicitly, leading to
the classical Laguerre orthogonal polynomials. In this case numerical values
to confirm the details of the spectrum can also be obtained by use of the
Liouville transformation; this leads to the Laguerre/Liouville Example 23, for
which the program is successful over a wide range of boundary conditions.

The Liouville transformation has also been applied to the Jacobi equation,
Example 16, to yield the Jacobi/Liouville Example 24.

The Liouville transformation is sometimes useful in other cases to put a
Sturm-Liouville differential equation into a form more suitable for numerical
computation; see in particular the Bessel Example 2.

\textbf{Parameters.} The reader is reminded that many of the examples involve
the choice of one or more parameters; the range of these parameters is given
when the numbered differential equation is displayed in xamples.x; if a choice
of parameter is made outside of the stated range the program may abort.

\section{Remarks on the individual examples.}

\begin{enumerate}
\item \textbf{Classical Legendre equation} (see \cite[Chapter IV]{T})
\[
-\left(  \left(  1-x^{2}\right)  y^{\prime}(x)\right)  ^{\prime}+\tfrac{1}%
{4}y(x)=\lambda y(x)\;\text{for all}\;x\in(-1,+1).
\]

Endpoint classification in $L^{2}(-1,+1)$:%
\[%
\begin{tabular}
[c]{cc}\hline
Endpoint & Classification\\\hline
$-1$ & LCNO\\
$+1$ & LCNO\\\hline
\end{tabular}
\]
$\quad$

For both endpoints the boundary condition functions $u,v$ are given by (note
that $u$ and $v$ are solutions of the Legendre equation for $\lambda=1/4$)%
\[
u(x)=1\quad\quad v(x)=\frac{1}{2}\ln\left(  \frac{1+x}{1-x}\right)
\;\text{for all}\;x\in(-1,+1).
\]

\begin{itemize}
\item[(i)] The Legendre polynomials are obtained by taking the principal
(Friedrichs) boundary condition at both endpoints $\pm1:$ enter
$A1=1,A2=0,\;B1=1,B2=0;$ \textit{i.e.} take the boundary condition function
$u$ at $\pm1$; eigenvalues: $\lambda_{n}=(n+1/2)^{2}\ ;$ $n=0,1,2,\cdots;$
eigenfunctions: Legendre polynomials $P_{n}(x)$.

\item[(ii)] Enter $A1=0,\;A2=1,\;B1=0,\;B2=1,$ \textit{i.e.} use the boundary
condition function $v$ at $\pm1$; eigenvalues: $\mu_{n};$ $n=0,1,2,\cdots$ but
no explicit formula is available; eigenfunctions are logarithmically unbounded
at $\pm1$.

\item[(iii)] Observe that $\mu_{n}<\lambda_{n}<\mu_{n+1};$ $n=0,1,2\cdots$.
\end{itemize}

\item \textbf{The Bessel equation }(see \cite[Chapter IV]{T})
\[
-y^{\prime\prime}(x)+\left(  \nu^{2}-1/4\right)  x^{-2}y(x)=\lambda
y(x)\;\text{for all}\;x\in(0,+\infty)
\]
with the parameter $\nu\in\lbrack0,+\infty).$ This is the Liouville form of
the classical Bessel equation.

Endpoint classification in $L^{2}(0,+\infty)$:%
\[%
\begin{tabular}
[c]{ccc}\hline
Endpoint & Parameter $\nu$ & Classification\\\hline
$0$ & For $\nu=1/2$ & R\\
$0$ & For all $\nu\in\lbrack0,1)$ but $\nu\neq1/2$ & LCNO\\
$0$ & For all $\nu\in\lbrack1,\infty)$ & LP\\\hline
$+\infty$ & For all $\nu\in\lbrack0,\infty)$ & LP\\\hline
\end{tabular}
\]

For endpoint $0$ and $\nu\in(0,1)$ but $\nu\neq1/2,$ the LCNO boundary
condition functions $u,v$ are determined by, for all$\;x\in(0,+\infty),$
\[%
\begin{tabular}
[c]{ccc}\hline
Parameter & $u$ & $v$\\\hline
$\nu\in(0,1)$ but $\nu\neq1/2$ & $x^{\nu+1/2}$ & $x^{-\nu+1/2}$\\
$\nu=0$ & $x^{1/2}$ & $x^{1/2}\ln(x)$\\\hline
\end{tabular}
\
\]

(a) Problems on $(0,1]$ with $y(1)=0$:

\noindent For $0\leq\nu<1,\nu\neq\frac{1}{2}:$ the Friedrichs case:
$A1=1,A2=0$ yields the classical Fourier-Bessel series; here $\lambda
_{n}=j_{\nu,n}^{2}$ where $\{j_{\nu,n}:n=0,1,2,\ldots\}$ are the zeros
(positive) of the Bessel function $J_{\nu}(\cdot).$

\noindent For $\nu\geq1;$ LP at $0$ so that there is a unique boundary value
problem with $\lambda_{n}=j_{\nu,n}^{2}$ as before.

(b) Problems on $[1,\infty)$ all have continuous spectrum on $[0,\infty)$:

\noindent For Dirichlet and Neumann boundary conditions there are no eigenvalues.

\noindent For $A1=A2=1$ at $1$ there is one isolated negative eigenvalue.

(c) Problems on $(0,\infty)$ all have continuous spectrum on $[0,\infty)$:

\noindent For $\nu\geq1$ there are no eigenvalues.

\noindent For $0\leq\nu<1$ the Friedrichs case is given by $A1=1,A2=0;$ there
are no eigenvalues.

\noindent For $\nu=0.45$ and $A1=10,A2=-1$ there is one isolated eigenvalue
near to the value $-175.57.$\noindent

\item \textbf{The Halvorsen equation}%
\[
-y^{\prime\prime}(x)=\lambda x^{-4}\exp(-2/x)y(x)\;\text{for all }%
x\in(0,+\infty)
\]

The endpoint classification in the weighted space $L^{2}((0,+\infty;x^{-4}%
\exp(-2/x))$:%
\[%
\begin{tabular}
[c]{cc}\hline
Endpoint & Classification\\\hline
$0$ & WR\\
$+\infty$ & LCNO
\end{tabular}
\]

For the endpoints $0$ and $+\infty$ in the WR and LCNO classification the
boundary condition functions $u,v$ are determined by%
\[%
\begin{tabular}
[c]{ccc}\hline
Endpoint & $u$ & $v$\\\hline
$0$ & $x$ & $1$\\
$+\infty$ & $1$ & $x$\\\hline
\end{tabular}
\]

Since this equation is WR at $0$ and LCNO at $+\infty$ the spectrum is
discrete and bounded below for all boundary conditions. However, this example
illustrates that even a R or WR endpoint can cause difficulties for
computation. The program fails on R at $0$; is successful for WR at $0$; is
successful for LCNO at $0.$

At $0$, the principal boundary condition entry is $A1=1,\ A2=0$; at $\infty$
with $u(x)=1,\ v(x)=x$ the principal boundary condition entry is also
$A1=1,\;A2=0,$ but note the interchange of the definitions of $u$ and $v$ at
these two endpoints.

\item \textbf{The Boyd equation}%
\[
-y^{\prime\prime}(x)-x^{-1}y(x)=\lambda y(x)\;\text{for all}\;x\in
(-\infty,0)\cup(0,+\infty).
\]

Endpoint classification in $L^{2}(-\infty,0)\cup$ $L^{2}(0,+\infty)$:%
\[%
\begin{tabular}
[c]{cc}\hline
Endpoint & Classification\\\hline
$-\infty$ & LP\\
$0-$ & LCNO\\
$0+$ & LCNO\\
$+\infty$ & LP\\\hline
\end{tabular}
\]

For both endpoints $0-$ and $0+$%
\[
u(x)=x\quad\quad v(x)=x\ln(\left|  x\right|  )\;\text{for all}\;x\in
(-\infty,0)\cup(0,+\infty).
\]

This equation arises in a model studying eddies in the atmosphere; see
\cite{B}. There is no explicit formula for the eigenvalues of any particular
boundary condition; eigenfunctions can be given in terms of Whittaker
functions; see \cite[Example 3]{BEZ}.

\item \textbf{The regularized Boyd equation}%
\[
-(p(x)y^{\prime}(x))^{\prime}+q(x)y(x)=\lambda w(x)y(x)\;\text{for all}%
\;x\in(-\infty,0)\cup(0,+\infty)
\]
where%
\[
p(x)=r(x)^{2}\quad q(x)=-r(x)^{2}\left(  \ln(\left|  x\right|  \right)
^{2}\quad w(x)=r(x)^{2}%
\]
with%
\[
r(x)=\exp\left(  -(x\ln(\left|  x\right|  )-x)\right)  \;\text{for all}%
\;x\in(-\infty,0)\cup(0,+\infty).
\]

Endpoint classification in $L^{2}(-\infty,0)\cup$ $L^{2}(0,+\infty)$:%
\[%
\begin{tabular}
[c]{cc}\hline
Endpoint & Classification\\\hline
$-\infty$ & LP\\
$0-$ & WR\\
$0+$ & WR\\
$+\infty$ & LP\\\hline
\end{tabular}
\]

This is a WR form of Example 4; the singularity at zero has been regularized
using quasi-derivatives. There is a close relationship between the examples 4
and 5; in particular they have the same eigenvalues - see \cite{AEZ}. For a
general discussion of regularization using non-principal solutions see
\cite{NZ}. For numerical results see \cite[Example 3]{BEZ}.

\item \textbf{The Sears-Titchmarsh equation}%
\[
-(xy^{\prime}(x))^{\prime}-xy(x)=\lambda x^{-1}y(x)\;\text{for all}%
\;x\in(0,+\infty).
\]

Endpoint classification in $L^{2}(0,\infty)$:%
\[%
\begin{tabular}
[c]{cc}\hline
Endpoint & Classification\\\hline
$0$ & LP\\
$+\infty$ & LCO\\\hline
\end{tabular}
\]

For the endpoint $+\infty$%
\[
u(x)=x^{-1/2}\left(  \cos(x)+\sin(x)\right)  \quad v(x)=x^{-1/2}\left(
\cos(x)-\sin(x)\right)  \;\text{for all}\;x\in(0,+\infty).
\]

This differential equation has one LP and one LCO endpoint. For details of
boundary value problems on $[1,\infty)$ see \cite[Example 4]{BEZ}. The
equation was studied originally in \cite[Chapter IV]{T}; but see \cite{ST}.

For problems on $[1,\infty)$ the spectrum is simple and discrete but unbounded
both above and below.

Numerical results are given in \cite[Example 4]{BEZ}.

\item \textbf{The BEZ equation}%
\[
-(xy^{\prime}(x))^{\prime}-x^{-1}y(x)=\lambda y(x)\;\text{for all}%
\;x\in(-\infty,0)\cup(0,+\infty).
\]

Endpoint classification in $L^{2}(-\infty,0)\cup$ $L^{2}(0,+\infty)$:%
\[%
\begin{tabular}
[c]{cc}\hline
Endpoint & Classification\\\hline
$-\infty$ & LP\\
$0-$ & LCO\\
$0+$ & LCO\\
$+\infty$ & LP\\\hline
\end{tabular}
\]

For both endpoints $0-$ and $0+$:%
\[
u(x)=\cos\left(  \ln(\left|  x\right|  )\right)  \quad\quad v(x)=\sin\left(
\ln(\left|  x\right|  )\right)  \;\text{for all}\;x\in(-\infty,0)\cup
(0,+\infty).
\]

This example is similar to the differential equation of Example 6. On the
interval $(0,1]$ there is a singularity at $0$ in LCO; the equation is R at 1.

For numerical results see \cite[Example 5]{BEZ}.

\item \textbf{The Laplace tidal wave equation}%
\[
-(x^{-1}y^{\prime}(x))^{\prime}+\left(  kx^{-2}+k^{2}x^{-1}\right)
y(x)=\lambda y(x)\;\text{for all}\;x\in(0,+\infty)
\]
where the parameter $k\in(-\infty,0)\cup(0,+\infty)$

Endpoint classification in $L^{2}(0,\infty)$:%
\[%
\begin{tabular}
[c]{cc}\hline
Endpoint & Classification\\\hline
$0$ & LCNO\\
$+\infty$ & LP\\\hline
\end{tabular}
\]

For the endpoint $0$:%
\[
u(x)=x^{2}\quad\quad v(x)=x-k^{-1}\;\text{for all}\;(0,+\infty).
\]

This equation is a particular case of the more general equation with this
name; for details and references see \cite{H}.

There are no representations for solutions of this differential equation in
terms of the well-known special functions. Thus to determine boundary
conditions at the LCNO endpoint $0$ use has to be made of maximal domain
functions; see the $u,\ v$ functions given above. Numerical results for some
boundary value problems and certain values of the parameter $k,$ are given in
\cite[Example 8]{BEZ}.

\item \textbf{The Latzko equation}%
\[
-((1-x^{7})y^{\prime}(x))^{\prime}=\lambda x^{7}y(x)\;\text{for all}%
\;x\in(0,1].
\]

Endpoint classification in $L^{2}(0,1]$:%
\[%
\begin{tabular}
[c]{cc}\hline
Endpoint & Classification\\\hline
$0$ & WR\\
$1$ & LCNO\\\hline
\end{tabular}
\]

For the endpoint $1$:%
\[
u(x)=1\quad\quad v(x)=-\ln(1-x)\;\text{for all}\;(0,1).
\]

This differential equation has a long and celebrated history; see \cite[Pages
43 to 45]{F}. There is a LCNO singularity at the endpoint $1$ which requires
the use of maximal domain functions; see the $u,\ v$ functions given above.
The endpoint $0$ is WR due to the fact that $w(0)=0$.

This example is similar in some respects to the Legendre equation of Example 1 above.

For numerical results see \cite[Example 7]{BEZ}.

\item \textbf{A weakly regular equation}%
\[
-(x^{1/2}y^{\prime}(x))^{\prime}=\lambda x^{-1/2}y(x)\;\text{for all}%
\;x\in(0,+\infty).
\]

Endpoint classification in $L^{2}(0,1]$:%
\[%
\begin{tabular}
[c]{cc}\hline
Endpoint & Classification\\\hline
$0$ & WR\\
$+\infty$ & LP\\\hline
\end{tabular}
\]

This is a devised example to illustrate the computational difficulties of
weakly regular problems.

The differential equation gives $p(0)=0$ and $w(0)=\infty$ but nevertheless
$0$ is a regular endpoint in the Lebesgue integral sense; however $0$ has to
be classified as weakly regular in the computational sense.

The Liouville normal form of this equation is the Fourier equation, see
Example 21 below; thus numerical results for this WR problem can be checked
against numerical results from (i) a R problem, (ii) the roots of
trigonometrical equations, and (iii) a LCNO problem (see below).

There are explicit solutions of this equation given by%
\[
\cos(2x^{1/2}\surd\lambda)\ \ ;\ \sin(2x^{1/2}\surd\lambda)/\surd\lambda.
\]

If $0$ is treated as a LCNO endpoint then $u,\ v$ boundary condition functions
are%
\[
u(x)=2x^{1/2}\quad\quad\ v(x)=1.
\]

The regular Dirichlet condition$\;y(0)=0$ is equivalent to the singular
condition $[y,u](0)=0$. Similarly the regular Neumann condition $(py^{\prime
})(0)=0$ is equivalent to the singular condition $[y,\ v](0)=0$.

The following indicated boundary value problems have the given explicit
formulae for the eigenvalues:
\begin{align*}
y(0)  &  =0\text{{ or }}[y,\ u](0)=0,\text{{ and }}y(1)=0\text{ gives}\\
\;\;\lambda_{n}  &  =((n+1)\pi)^{2}/4\ (n=0,1,...)
\end{align*}%
\begin{align*}
(py^{\prime})(0)  &  =0\text{{ or }}[y,v](0)=0,\text{{ and }}(py^{\prime
})(1)=0\text{ gives}\\
\;\;\lambda_{n}  &  =\left(  (n+\tfrac{1}{2})\pi\right)  ^{2}/4\ (n=0,1,...).
\end{align*}

\item \textbf{The Plum equation}%
\[
-(y^{\prime}(x))^{\prime}+100\cos^{2}(x)y(x)=\lambda y(x)\;\text{for
all}\;x\in(-\infty,+\infty).
\]

Endpoint classification in $L^{2}(-\infty,+\infty)$:%
\[%
\begin{tabular}
[c]{cc}\hline
Endpoint & Classification\\\hline
$-\infty$ & LP\\
$+\infty$ & LP\\\hline
\end{tabular}
\]

Plum \cite{P} computed the first seven eigenvalues for periodic eigenvalues on
the interval $[0,\pi],$\textit{i.e.}%
\[
y(0)=y(\pi)\quad\quad\quad y^{\prime}(0)=y^{\prime}(\pi),
\]
using a numerical homotopy method together with interval arithmetic, and
obtained rigorous bounds for these seven computed eigenvalues. In double
precision the SLEIGN2 computed eigenvalues are in good agreement with these
guaranteed bounds.

\item \textbf{The Mathieu equation}%
\[
-y^{\prime\prime}(x)+2k\cos(2x)y(x)=\lambda y(x)\;\text{for all}\;x\in
(-\infty,+\infty)
\]
where that parameter $k\in(-\infty,0)\cup(0,+\infty).$

Endpoint classification in $L^{2}(-\infty,+\infty)$:%
\[%
\begin{tabular}
[c]{cc}\hline
Endpoint & Classification\\\hline
$-\infty$ & LP\\
$+\infty$ & LP\\\hline
\end{tabular}
\]

The classical Mathieu equation has a celebrated history and voluminous
literature. There are no eigenvalues for this problem on $(-\infty,+\infty)$.
There may be one negative eigenvalue of the problem on $[0,\infty)$ depending
on the boundary condition at the endpoint $0$. The continuous (essential)
spectrum is the same for the whole line or half-line problems and consists of
an infinite number of disjoint closed intervals. The endpoints of these - and
thus the spectrum of the problem - can be characterized in terms of periodic
and semi-periodic eigenvalues of Sturm-Liouville problems on the compact
interval $[0,2\pi]$; these can be computed with SLEIGN2.

The above remarks also apply to the general Sturm-Liouville equation with
periodic coefficients of the same period; the so-called Hill's equation.

Of special interest is the starting point of the continuous spectrum - this is
also the oscillation number of the equation. For the Mathieu equation
($p=1,q=\cos(x),w=1$) on both the whole line and the half line it is
approximately -0.378; this result may be obtained by computing the first
eigenvalue $\lambda_{0}$ of the periodic problem on the interval $[0,2\pi]$.

\item \textbf{The hydrogen atom equation}

It is convenient to take this equation in two forms:%
\begin{equation}
-y^{\prime\prime}(x)+(kx^{-1}+hx^{-2})y(x)=\lambda y(x)\;\text{for all}%
\;x\in(0,+\infty) \tag{1}%
\end{equation}
where the two independent parameters $h\in\lbrack-1/4,+\infty)$ and
$k\in\mathbb{R},$ and%
\begin{equation}
-y^{\prime\prime}(x)+(kx^{-1}+hx^{-2}+1)y(x)=\lambda y(x)\;\text{for
all}\;x\in(0,+\infty) \tag{2}%
\end{equation}
where the two independent parameters $h\in(-\infty,-1/4)$ and $k\in\mathbb{R}.$

Note that form (2) is introduced as a device to aid the numerical computations
in the difficult LCO case; it forces the boundary value problem to have a
non-negative eigenvalue.

Endpoint classification, for both forms (1) and (2), in $L^{2}(0,+\infty)$:%
\[%
\begin{tabular}
[c]{cccc}\hline
Endpoint & Form & Parameters & Classification\\\hline
$0$ & 1 & $h=k=0$ & R\\
$0$ & 1 & $h=0,k\in\mathbb{R}\,\backslash\,\{0\}$ & LCNO\\
$0$ & 1 & $-1/4\leq h<3/4,h\neq0,k\in\mathbb{R}$ & LCNO\\
$0$ & 1 & $h\geq3/4,k\in\mathbb{R}$ & LP\\
$0$ & 2 & $h<-1/4,k\in\mathbb{R}$ & LCO\\\hline
$+\infty$ & 1 and 2 & $h,k\in\mathbb{R}$ & LP\\\hline
\end{tabular}
\]

This is the two parameter version of the classical one-dimensional equation
for quantum modelling of the hydrogen atom; see \cite[Section 10]{J}.

For form (1) and all $h,k$ there are no positive eigenvalues; form (2) is best
considered in the single LCO case when some eigenvalues are positive; in form
(1) there is a continuous spectrum on $[0,\infty)$; in form (2) there is a
continuous spectrum on $[1,\infty).$

If $k=0$ the equation reduces to Bessel, see Example 2 above with $h=\nu
^{2}-1/4$.

\noindent\textbf{Results for form (1)}

In all cases below $\rho$ is defined by
\[
\rho:=(h+1/4)^{1/2}\text{ for }h\geq-1/4.
\]

\begin{enumerate}
\item For $h\geq3/4$ and $k\geq0$ no boundary conditions are required; there
is at most one negative eigenvalue and $\lambda=0$ may be an eigenvalue; for
$h\geq3/4$ and $k<0$ there are infinitely many negative eigenvalues given by
\[
\lambda_{n}=\frac{-k^{2}}{(2n+2\rho+1)^{2}},\ \rho=(h+1/4)^{1/2}%
>0,\ n=0,1,2,3,\ldots
\]
and $\lambda=0$ is not an eigenvalue.

\item For $h=0$ and $k\in\mathbb{R}\,\backslash\,\{0\}$ a boundary condition
is required at $0$ for which
\[
u(x)=x\quad\quad\ v(x)=1+k\,x\ln(x).
\]
For some computed eigenvalues see \cite{BEZ} and \cite[Section 10]{J}.

\item For $-1/4<h<3/4,$ \textit{i.e.} $0<\rho<1,$ and $h\neq0,$ \textit{i.e.}
$\rho\not =1/2$, then a boundary condition is required at $0$ for which, for
all $x\in(0,+\infty),$%
\[
u(x)=x^{\frac{1}{2}+\rho}\quad\quad v(x)=x^{\frac{1}{2}-\rho}+\frac{k}%
{1-2\rho}x^{\frac{3}{2}-\rho};
\]

the following results hold for the non-Friedrichs boundary condition
$[y,v](0)=0$, \textit{i.e.} $A1=0,A2=1$:

\begin{itemize}
\item[(i)] $k>0,\ 0<\rho<1/2$ there are no negative eigenvalues

\item[(ii)] $k>0,\ 1/2<\rho<1$ there is exactly one negative eigenvalue given
by
\[
\lambda_{0}=\frac{-k^{2}}{(2\rho-1)^{2}}%
\]

\item[(iii)] if $k<0,\ 0<\rho<1/2$ there are infinitely many negative
eigenvalues given by
\[
\lambda_{n}=\frac{-k^{2}}{(2n-2\rho+1)^{2}},\ n=0,1,2,3,\ldots
\]

\item[(iv)] if $k<0,\ 1/2<\rho<1$ there are infinitely many negative
eigenvalues given by
\[
\lambda_{n}=\frac{-k^{2}}{(2n-2\rho+3)^{2}}\,,\ n=0,1,2,3,\ldots
\]

\item[(v)] for $k=0$ and $(A1)(A2)<0$ there is exactly one negative eigenvalue
given by:
\[
\lambda_{0}=\frac{4\left(  A1\right)  \Gamma(1+\rho)}{\left(  A2\right)
\Gamma(1-\rho)^{1/\rho}}.
\]
\end{itemize}

\item For $h=-1/4,\ k\in R$ , the LCNO classification at $0$ prevails and a
boundary condition is required for which, for all $x\in(0,+\infty),$%
\[
u(x)=x^{1/2}+kx^{3/2}\quad\quad v(x)=2x^{1/2}+\left(  x^{1/2}+kx^{3/2}\right)
\ln(x).
\]
\noindent For $k=0$ and $(A1)(A2)<0$ there is exactly one negative eigenvalue
given by:
\[
\lambda_{0}=-c\exp(2(A1)/A2),\quad c=4\exp(4-2\gamma)
\]
where $\gamma$ is Euler's constant: $\gamma=0.5772156649\ldots$.

\noindent\noindent\textbf{Results for form (2)}

\item For $h<-1/4,\ k\in R$, the equation is LCO at $0$ (recall that we added
$1$ to the coefficient $q(\cdot)$ for this case, thus moving the start of the
continuous spectrum from $0$ to $1$) for which, defining%
\[
\sigma:=(-h-1/4)^{1/2},
\]
then, for all $x\in(0,+\infty),$%
\begin{align*}
u(x)  &  =x^{1/2}\left[  (1-(4h)^{-1}kx)\cos(\sigma\ln(x))+k\sigma
x\sin(\sigma\ln(x))/2\right] \\
v(x)  &  =x^{1/2}\left[  (1-(4h)^{-1}kx)\sin(\sigma\ln(x))+k\sigma
x\cos(\sigma\ln(x))/2\right]  ;
\end{align*}

\begin{itemize}
\item[(i)] when $k=0$ this equation reduces to the Krall equation Example 20
below (but note that the notation is different).

\item[(ii)] When $k\not =0$ explicit formulas for the eigenvalues are not
available; however we report here on the qualitative properties of the
spectrum for any boundary condition at $0$:

$(\alpha)$ for all $k\in R$ there are infinitely many negative eigenvalues
tending exponentially to $-\infty$

$(\beta)$ for $k>0$ there are only a finite number of eigenvalues in any
bounded interval, in particular they do not accumulate at $1$

$(\gamma)$ for $k\leq0$ the eigenvalues accumulate also at $1$.

$(\delta)$ for $k=0$ and $(A1)(A2)<0$ there is exactly one negative eigenvalue
given by:
\[
\lambda_{0}=\frac{4\left(  A1\right)  \Gamma(1+\rho)}{\left(  A2\right)
\Gamma(1-\rho)^{1/\rho}}.
\]
\end{itemize}

Most of these results are due to J\"{o}rgens, see \cite[Section 10]{J}; a few
new results were established by the authors.
\end{enumerate}

\item \textbf{The Marletta equation}%
\[
-y^{\prime\prime}(x)+\frac{3(x-31)}{4(x+1)(x+4)^{2}}y(x)=\lambda
y(x)\;\text{for all}\;x\in\lbrack0,+\infty).
\]

Endpoint classification in $L^{2}(0,+\infty)$:%
\[%
\begin{tabular}
[c]{cc}\hline
Endpoint & Classification\\\hline
$0$ & R\\
$+\infty$ & LP\\\hline
\end{tabular}
\]
Since $q(x)\rightarrow0$ as $x\rightarrow\infty$ the continuous spectrum
consists of $[0,\infty)$ and every negative number is an eigenvalue for some
boundary condition at $0.$

For the boundary condition $A1=5,A2=8$ there is a negative eigenvalue
$\lambda_{0}$ near $-1.185.$ However the equation with $\lambda=0$ has a
solution
\[
y(x)=\frac{1-x^{2}}{(1+x/4)^{5/2}}\text{ for all }x\in\lbrack0,\infty)
\]
that satisfies this boundary condition which is NOT in $L^{2}(0,\infty)$ but
is ``nearly'' in this space. This solution deceives SLEIGN and SLEIGN2 in
single precision, and SLEDGE in double precision, into reporting $\lambda=0$
as a second eigenvalue; in double precision SLEIGN and SLEIGN2 correctly
report that $\lambda_{0}$ is the only eigenvalue, and SLEIGN2 reports the
start of the continuous spectrum at $0.$

Additional details of this example are to be found in the Marletta
certification report on SLEIGN (not SLEIGN2) \cite{M}.

\item \textbf{The harmonic oscillator equation}%
\[
-y^{\prime\prime}(x)+x^{2}y(x)=\lambda y(x)\;\text{for all}\;x\in
(-\infty,+\infty).
\]

Endpoint classification in $L^{2}(-\infty,+\infty)$:%
\[%
\begin{tabular}
[c]{cc}\hline
Endpoint & Classification\\\hline
$-\infty$ & LP\\
$+\infty$ & LP\\\hline
\end{tabular}
\
\]

This is another classical equation; it is also the Liouville normal form of
the differential equation for the Hermite orthogonal polynomials. On the whole
real line the boundary value problem requires no boundary conditions at the
endpoints of $\pm\infty$. Thus there is a unique self-adjoint extension with
discrete spectrum given by :
\[
\{\lambda_{n}=2n+1;\;n=0,1,2,...\}.
\]
For a classical treatment see \cite[Chapter IV, Section 2]{T}.

\item \textbf{The Jacobi equation}%
\[
-\left(  (1-x)^{\alpha+1}(1+x)^{\beta+1}y^{\prime}(x)\right)  ^{\prime
}=\lambda(1-x)^{\alpha}(1+x)^{\beta}y(x)\;\text{for all}\;x\in(-1,+1)
\]
where the parameters $\alpha,\beta\in(-\infty,+\infty).$

Endpoint classification in the weighted space $L^{2}((-1,+1);(1-x)^{\alpha
}(1+x)^{\beta}))$:%
\[%
\begin{tabular}
[c]{ccc}\hline
Endpoint & Parameter & Classification\\\hline
$-1$ & $\beta\leq-1$ & LP\\
$-1$ & $-1<\beta<0$ & WR\\
$-1$ & $0\leq\beta<1$ & LCNO\\
$-1$ & $1\leq\beta$ & LP\\\hline
\end{tabular}
\quad\quad%
\begin{tabular}
[c]{ccc}\hline
Endpoint & Parameter & Classification\\\hline
$+1$ & $\alpha\leq-1$ & LP\\
$+1$ & $-1<\alpha<0$ & WR\\
$+1$ & $0\leq\alpha<1$ & LCNO\\
$+1$ & $1\leq\alpha$ & LP\\\hline
\end{tabular}
\]

For the endpoint $-1$ and for the WR and LCNO cases the boundary condition
functions $u,v$ are determined by%
\[%
\begin{tabular}
[c]{ccc}\hline
Parameter & $u$ & $v$\\\hline
$-1<\beta<0$ & $(1+x)^{-\beta}$ & $1$\\
$\beta=0$ & $1$ & $\ln\left(  \dfrac{1+x}{1-x}\right)  $\\
$0<\beta<1$ & $1$ & $(1+x)^{-\beta}$\\\hline
\end{tabular}
.
\]

For the endpoint $+1$ and for the WR and LCNO cases the boundary condition
functions $u,v$ are determined by%
\[%
\begin{tabular}
[c]{ccc}\hline
Parameter & $u$ & $v$\\\hline
$-1<\alpha<0$ & $(1-x)^{-\alpha}$ & $1$\\
$\alpha=0$ & $1$ & $\ln\left(  \dfrac{1+x}{1-x}\right)  $\\
$0<\alpha<1$ & $1$ & $(1-x)^{-\alpha}$\\\hline
\end{tabular}
.
\]

To obtain the classical Jacobi orthogonal polynomials it is necessary to take
$-1<\alpha,\,\beta$; then note the required boundary conditions:

Endpoint $-1$:%
\[%
\begin{tabular}
[c]{cc}\hline
Parameter & Boundary condition\\\hline
$-1<\beta<0$ & $(py^{\prime})(-1)=0\;$or$\;[y,v](-1)=0$\\
$0\leq\beta<1$ & $[y,u](-1)=0$\\\hline
\end{tabular}
\]

Endpoint $+1$:%
\[%
\begin{tabular}
[c]{cc}\hline
Parameter & Boundary condition\\\hline
$-1<\alpha<0$ & $(py^{\prime})(+1)=0\;$or$\;[y,v](+1)=0$\\
$0\leq\alpha<1$ & $[y,u](+1)=0$\\\hline
\end{tabular}
\]
\newline For the classical Jacobi orthogonal polynomials the eigenvalues are
given by:
\[
\lambda_{n}=n(n+\alpha+\beta+1)\;\text{for}\;n=0,1,2,\ldots
\]
and this explicit formula can be used to give an independent check on the
accuracy of the results from the SLEIGN2 code.

It is interesting to note that the required boundary condition for these
Jacobi polynomials is the Friedrichs condition in the LCNO case but not in the
WR case.

\item \textbf{The rotation Morse oscillator equation}%
\[
-y^{\prime\prime}(x)+(2x^{-2}-2000(2e(x)-e(x)^{2}))y(x)=\lambda
y(x)\;\text{for all}\;x\in(0,+\infty)
\]
where%
\[
e(x)=\exp(-1.7(x-1.3))\;\text{for all}\;x\in(0,+\infty).
\]

Endpoint classification in the space $L^{2}(0,+\infty)$%
\[%
\begin{tabular}
[c]{cc}\hline
Endpoint & Classification\\\hline
$0$ & LP\\
$+\infty$ & LP\\\hline
\end{tabular}
\
\]

This classical problem has continuous spectrum on $[0,\infty)$ and exactly 26
negative eigenvalues. Enter NUMEIG1 = 0, NUMEIG2 = 28 and observe the 26
eigenvalues and the start of the continuous spectrum at 0.

\item \textbf{The Dunsch equation}%
\[
-\left(  (1-x^{2})y^{\prime}(x)\right)  ^{\prime}+\left(  \frac{2\alpha^{2}%
}{(1+x)}+\frac{2\beta^{2}}{(1-x)}\right)  y(x)=\lambda y(x)\;\text{for
all}\;x\in(-1,+1)
\]
where the independent parameters $\alpha,\beta\in\lbrack0,+\infty).$

Boundary value problems for this differential equation are discussed in
\cite[Chapter VIII, Pages 1515-20]{DS}.

Endpoint classification in the space $L^{2}(-1,+1)$ for $-1$:%
\[%
\begin{tabular}
[c]{cc}\hline
Parameter & Classification\\\hline
$0\leq\alpha<1/2$ & LCNO\\
$1/2\leq\alpha$ & LP\\\hline
\end{tabular}
\]

Endpoint classification in the space $L^{2}(-1,+1)$ for $+1$:%
\[%
\begin{tabular}
[c]{cc}\hline
Parameter & Classification\\\hline
$0\leq\beta<1/2$ & LCNO\\
$1/2\leq\beta$ & LP\\\hline
\end{tabular}
\]

For the LCNO cases the boundary condition functions $u,v$ are given by%
\[%
\begin{tabular}
[c]{cccc}\hline
Endpoint & Parameter & $u$ & $v$\\\hline
$-1$ & $\alpha=0$ & $1.0$ & $\dfrac{1}{2}\ln\left(  \dfrac{1+x}{1-x}\right)
$\\
$-1$ & $0<\alpha<1/2$ & $(1+x)^{\alpha}$ & $(1+x)^{-\alpha}$\\
$+1$ & $\beta=0$ & $1.0$ & $\dfrac{1}{2}\ln\left(  \dfrac{1+x}{1-x}\right)
$\\
$+1$ & $0<\beta<1/2$ & $(1-x)^{\beta}$ & $(1-x)^{-\beta}$\\\hline
\end{tabular}
\]

Note that these $u$ and $v$ are not solutions of the differential equation but
maximal domain functions. In \cite[Page 1519]{DS} it is stated that the
boundary value problem determined by the boundary conditions
\[
\lbrack y,u](-1)=0=[y,u](1)
\]
has eigenvalues given by the explicit formula
\[
\lambda_{n}=(n+\alpha+\beta+1)(n+\alpha+\beta)\;\text{for}\;n=0,1,2,\ldots
\]

\item \textbf{The Donsch equation}%
\[
-\left(  (1-x^{2})y^{\prime}(x)\right)  ^{\prime}+\left(  \frac{-2\gamma^{2}%
}{(1+x)}+\frac{2\beta^{2}}{(1-x)}\right)  y(x)=\lambda y(x)\;\text{for
all}\;x\in(-1,+1)
\]
where the independent parameters $\gamma,\beta\in\lbrack0,+\infty).$

Endpoint classification in the space $L^{2}(-1,+1)$ for $-1$:%
\[%
\begin{tabular}
[c]{cc}\hline
Parameter & Classification\\\hline
$\gamma=0$ & LCNO\\
$0<\gamma$ & LCO\\\hline
\end{tabular}
\]

Endpoint classification in the space $L^{2}(-1,+1)$ for $+1$:%
\[%
\begin{tabular}
[c]{cc}\hline
Parameter & Classification\\\hline
$0\leq\beta<1/2$ & LCNO\\
$1/2\leq\beta$ & LP\\\hline
\end{tabular}
\]

For these LCNO/LCO cases the boundary condition functions $u,v$ are given by%
\[%
\begin{tabular}
[c]{cccc}\hline
Endpoint & Parameter & $u$ & $v$\\\hline
$-1$ & $\gamma=0$ & $1$ & $\dfrac{1}{2}\ln\left(  \dfrac{1+x}{1-x}\right)  $\\
$-1$ & $0<\gamma$ & $\cos(\gamma\ln(1+x))$ & $\sin(\gamma\ln(1+x))$\\
$+1$ & $\beta=0$ & $1$ & $\dfrac{1}{2}\ln\left(  \dfrac{1+x}{1-x}\right)  $\\
$+1$ & $0<\beta<1/2$ & $(1-x)^{\beta}$ & $(1-x)^{-\beta}$\\\hline
\end{tabular}
\
\]

This is a modification of Example 18 above which illustrates an LCNO/LCO mix
obtained by replacing $\alpha$ with $i\gamma$; this changes the singularity at
$-1$ from LCNO to LCO.

Again these $u$ and $v$ are not solutions of the differential equation but
maximal domain functions.

\item \textbf{The Krall equation}%
\[
-y^{\prime\prime}(x)+(1-(k^{2}+1/4)x^{-2})y(x)=\lambda y(x)\;\text{for
all}\;x\in(0,+\infty)
\]
where the parameter $k\in(0,+\infty).$

Endpoint classification in the space $L^{2}(0,+\infty)$:%
\[%
\begin{tabular}
[c]{cc}\hline
Endpoint & Classification\\\hline
$0$ & LCO\\
$+\infty$ & LP\\\hline
\end{tabular}
\]

This example should be seen as a special case of the Bessel Example 2 above;
solutions can be obtained in terms of the modified Bessel functions.

To help with the computations for this example the spectrum is translated by a
term $+1$; this simple device is used for numerical convenience.

For problems with separated boundary conditions at endpoints $0$ and $\infty$
there is a continuous spectrum on $[1,\infty)$ with a discrete (and simple)
spectrum on $(-\infty,1)$. This discrete spectrum has cluster points at both
$-\infty$ and $1$.

For the LCO endpoint at $0$ the boundary condition functions are given by%
\[
u(x)=x^{1/2}\cos(k\ln(x))\quad\quad v(x)=x^{1/2}\sin(k\ln(x)).
\]

For the boundary value problem with boundary condition $[y,u](0)=0$ the
eigenvalues are given explicitly by:

(i) suppose $\Gamma(1+i)=\alpha+i\beta$ and $\mu>0$ satisfies $\tan\left(
\ln(\frac{1}{2}\mu)\right)  =-\alpha/\beta$

(ii) $\theta=\operatorname{Im}(\log(\Gamma(1+i)))$

(iii) $\ln(\frac{1}{2}\mu)=\tfrac{1}{2}\pi+\theta+s\pi\;$for$\;\;s=0,\pm1,\pm2,...$

(iv) $\mu_{s}^{2}=\left(  2\exp(\theta+\frac{1}{2}\pi)\right)  ^{2}\exp
(2s\pi)\;\,s=0,\pm1,\pm2,...$

\noindent then the eigenvalues are $\lambda_{n}=-\mu_{-(n+1)}^{2}%
+1\ (n=0,\pm1,\pm2,...)$.

SLEIGN2 can compute only six of these eigenvalues in a normal UNIX server,
even in double precision, $\lambda_{-3}$ to $\lambda_{2}$; other eigenvalues
are, numerically, too close to 1 or too close to $-\infty$. Here we list these
SLEIGN2 computed eigenvalues in double precision in a normal UNIX server and
compare them with the same eigenvalues computed from the transcendental
equation; for the problem on $(0,\infty)$ with $k=1$ and $A1=1.0,\ A2=0.0$.
\[%
\begin{array}
[c]{cccc}%
\text{NUMEIG} & \text{eig from SLEIGN2} & \text{eig from trans. equ.} &
\text{iflag}\\
-3 & -276,562.5 & -14,519,130 & 4\\
-2 & -27,114.48 & -27,114.67 & 2\\
-1 & -49.62697 & -49.63318 & 2\\
0 & 0.9054452 & 0.9054454 & 1\\
1 & 0.9998234 & 0.9998234 & 1\\
2 & 0.9999997 & 0.9999997 & 3
\end{array}
.
\]

\item \textbf{The Fourier equation}%
\[
-y^{\prime\prime}(x)=\lambda y(x)\;\text{for all}\;x\in(-\infty,+\infty)
\]

Endpoint classification in $L^{2}(-\infty,+\infty)$:%
\[%
\begin{tabular}
[c]{cc}\hline
Endpoint & Classification\\\hline
$-\infty$ & LP\\
$+\infty$ & LP\\\hline
\end{tabular}
\]

This is a simple constant coefficient equation whose eigenvalues, for any
self-adjoint boundary condition, can be characterized in terms of a
transcendental equation involving only trigonometric functions.

\item \textbf{The Laguerre equation}%
\[
-(x^{\alpha+1}\exp(-x)y^{\prime}(x))^{\prime}=\lambda x^{\alpha}%
\exp(-x)y(x)\;\text{for all}\;x\in(0,+\infty)
\]
where the parameter $\alpha\in(-\infty,+\infty).$

Endpoint classification in the weighted space $L^{2}((0,+\infty);x^{\alpha
}\exp(-x))$:%
\[%
\begin{tabular}
[c]{ccc}\hline
Endpoint & Parameter & Classification\\\hline
$0$ & $\alpha\leq-1$ & LP\\
$0$ & $-1<\alpha<0$ & WR\\
$0$ & $0\leq\alpha<1$ & LCNO\\
$0$ & $1\leq\alpha$ & LP\\
$+\infty$ & $\alpha\in(-\infty,+\infty)$ & LP\\\hline
\end{tabular}
\]

For these WR/LCNO cases the boundary condition functions $u,v$ are given by:%
\[%
\begin{tabular}
[c]{cccc}\hline
Endpoint & Parameter & $u$ & $v$\\\hline
$0$ & $-1<\alpha<0$ & $x^{-\alpha}$ & $1$\\
$0$ & $\alpha=0$ & $1$ & $\ln(x)$\\
$0$ & $0<\alpha<1$ & $1$ & $x^{-\alpha}$\\\hline
\end{tabular}
\]

This is the classical form of the differential equation which for parameter
$\alpha>-1$ produces the classical Laguerre polynomials as eigenfunctions; for
the boundary condition $[y,1](0)=0$ at $0$, when required, the eigenvalues are
then (remarkably!) independent of $\alpha$ and given by $\lambda
_{n}=n\;(n=0,1,2,...)$; see \cite[Chapter 22, Section 22.6]{AB}.

SLEIGN2 does not compute eigenvalues well with this differential equation on
$(0,\infty)$, with the code in a UNIX server; this appears to be due to
numerical problems resulting from the exponentially small coefficients;
however, see Example 23 below.

\item \textbf{The Laguerre/Liouville equation}%
\[
-y^{\prime\prime}(x)+\left(  \frac{\alpha^{2}-1/4}{x^{2}}-\frac{\alpha+1}%
{2}+\frac{x^{2}}{16}\right)  y(x)=\lambda y(x)\;\text{for all}\;x\in
(0,+\infty)
\]
where the parameter $\alpha\in(-\infty,+\infty).$

Endpoint classification in the space $L^{2}(0,+\infty)$:%
\[%
\begin{tabular}
[c]{ccc}\hline
Endpoint & Parameter & Classification\\\hline
$0$ & $\alpha\leq-1$ & LP\\
$0$ & $-1<\alpha<1,$ but $\alpha^{2}\neq1/4$ & LCNO\\
$0$ & $\alpha^{2}=1/4$ & R\\
$0$ & $1\leq\alpha$ & LP\\
$+\infty$ & $\alpha\in(-\infty,+\infty)$ & LP\\\hline
\end{tabular}
\]

For these WR/LCNO cases the boundary condition functions $u,v$ are given by:%
\[%
\begin{tabular}
[c]{cccc}\hline
Endpoint & Parameter & $u$ & $v$\\\hline
$0$ & $-1<\alpha<0$ but $\alpha\neq-1/2$ & $x^{\frac{1}{2}-\alpha}$ &
$x^{\frac{1}{2}+\alpha}$\\
$0$ & $\alpha=-1/2$ & $x$ & $1$\\
$0$ & $\alpha=0$ & $x^{1/2}$ & $x^{1/2}\ln(x)$\\
$0$ & $0<\alpha<1$ but $\alpha\neq1/2$ & $x^{\frac{1}{2}+\alpha}$ &
$x^{\frac{1}{2}-\alpha}$\\
$0$ & $\alpha=1/2$ & $x$ & $1$\\\hline
\end{tabular}
\]

This is the Liouville normal form of the Laguerre equation; the two forms are
unitarily equivalent so that the spectrum and the eigenfunctions of equivalent
boundary value problems are identical. This Liouville form is more suitable
for eigenvalue computations in contrast to the previous example.

The Laguerre polynomials are produced as eigenfunctions only when $\alpha>-1$.
For $\alpha\geq1$ the LP condition holds at $0$. For $0\leq\alpha<1$ the
appropriate boundary condition is the Friedrichs condition: $[y,u](0)=0;$ for
$-1<\alpha<0$ use the non-Friedrichs condition: $[y,\ v](0)=0$. In all these
cases $\lambda_{n}=n$ for $n=0,1,2,...$.

\item \textbf{The Jacobi/Liouville equation}%
\[
-y^{\prime\prime}(x)+q(x)y(x)=\lambda y(x)\;\text{for all}\;x\in(-\pi
/2,+\pi/2)
\]
where the coefficient $q$ is given by, for all$\;x\in(-\pi/2,+\pi/2),$%
\[
q(x)=\frac{\beta^{2}-1/4}{4\tan^{2}((x+\pi)/2)}+\frac{\alpha^{2}-1/4}%
{4\tan^{2}((x-\pi)/2)}-\frac{4\alpha\beta+4\beta+4\alpha+3}{8}.
\]
Here the parameters $\alpha,\beta\in(-\infty+,\infty).$

Endpoint classification in the space $L^{2}(-\pi/2,+\pi/2)$:%
\[%
\begin{tabular}
[c]{ccc}\hline
Endpoint & Parameter & Classification\\\hline
$-\pi/2$ & $\beta\leq-1$ & LP\\
$-\pi/2$ & $-1<\beta<1$ but $\beta^{2}\neq1/4$ & LCNO\\
$-\pi/2$ & $\beta^{2}=1/4$ & R\\
$-\pi/2$ & $1\leq\beta$ & LP\\\hline
\end{tabular}
\]%
\[%
\begin{tabular}
[c]{ccc}\hline
Endpoint & Parameter & Classification\\\hline
$+\pi/2$ & $\alpha\leq-1$ & LP\\
$+\pi/2$ & $-1<\alpha<1$ but $\alpha^{2}\neq1/4$ & LCNO\\
$+\pi/2$ & $\alpha^{2}=1/4$ & R\\
$+\pi/2$ & $1\leq\alpha$ & LP\\\hline
\end{tabular}
\]

For the endpoint $-\pi/2$ and for LCNO cases the boundary condition functions
$u,v$ are determined by, here $b(x)=2\tan^{-1}(1)+x$ for all $x\in(-\pi
/2,+\pi/2),$%
\[%
\begin{tabular}
[c]{ccc}\hline
Parameter & $u$ & $v$\\\hline
$-1<\beta<0$ & $b(x)^{\frac{1}{2}-\beta}$ & $b(x)^{\frac{1}{2}+\beta}$\\
$\beta=0$ & $\sqrt{b(x)}$ & $\sqrt{b(x)}\ln(b(x))$\\
$0<\beta<1$ & $b(x)^{\frac{1}{2}+\beta}$ & $b(x)^{\frac{1}{2}-\beta}$\\\hline
\end{tabular}
\]

For the endpoint $+\pi/2$ and for LCNO cases the boundary condition functions
$u,v$ are determined by, here $a(x)=2\tan^{-1}(1)-x$ for all $x\in(-\pi
/2,+\pi/2),$%
\[%
\begin{tabular}
[c]{ccc}\hline
Parameter & $u$ & $v$\\\hline
$-1<\alpha<0$ & $a(x)^{\frac{1}{2}-\alpha}$ & $a(x)^{\frac{1}{2}+\alpha}$\\
$\alpha=0$ & $\sqrt{a(x)}$ & $\sqrt{a(x)}\ln(a(x))$\\
$0<\alpha<1$ & $a(x)^{\frac{1}{2}+\alpha}$ & $a(x)^{\frac{1}{2}-\alpha}%
$\\\hline
\end{tabular}
\]

This is the Liouville normal form of the Jacobi equation of Example 16.

The classical Jacobi orthogonal polynomials are produced only when both
$\alpha,\beta>-1.$ For $\alpha,\beta>+1$ the LP condition holds and no
boundary condition is required to give the polynomials. If $-1<\alpha,\beta<1$
then the LCNO condition holds and boundary conditions are required to produce
the Jacobi polynomials; these conditions are as follows:

Endpoint $-\pi/2$%
\[%
\begin{tabular}
[c]{cc}\hline
Parameter & Boundary condition\\\hline
$-1<\beta<0$ & $[y,v](-\pi/2)=0$\\
$0\leq\beta<1$ & $[y,u](-\pi/2)=0$\\\hline
\end{tabular}
\]

Endpoint $+\pi/2$%
\[%
\begin{tabular}
[c]{cc}\hline
Parameter & Boundary condition\\\hline
$-1<\alpha<0$ & $[y,v](+\pi/2)=0$\\
$0\leq\alpha<1$ & $[y,u](+\pi/2)=0$\\\hline
\end{tabular}
\]

Recall from Example 16 for the classical orthogonal Jacobi polynomials the
eigenvalues are given explicitly by:%
\[
\lambda_{n}=n(n+\alpha+\beta+1)\;\text{for}\;n=0,1,2,\ldots
\]

\item \textbf{The Meissner equation}%
\[
-y^{\prime\prime}(x)=\lambda w(x)y(x)\;\text{for all}\;x\in(-\infty,+\infty)
\]
where the weight coefficient $w$ is defined by%
\begin{align*}
w(x)  &  =1\;\text{for all}\;x\in(-\infty,0]\\
&  =9\;\text{for all}\;x\in(0,+\infty).
\end{align*}

Endpoint classification in the space $L^{2}(-\infty,+\infty)$:%
\[%
\begin{tabular}
[c]{cc}\hline
Endpoint & Classification\\\hline
$-\infty$ & LP\\
$+\infty$ & LP\\\hline
\end{tabular}
\]

This equation arose in a model of a one dimensional crystal. For this constant
coefficient equation with a weight function which has a jump discontinuity the
eigenvalues can be characterized as roots of a transcendental equation
involving only trigonometrical and inverse trigonometrical functions. There
are infinitely many simple eigenvalues and infinitely many double ones for the
periodic case; they are given by:

\textbf{Periodic boundary conditions on} $(-1/2,+1/2)$\textbf{, }\textit{i.e.}%
\[
y(-1/2)=y(+1/2)\quad\quad y^{\prime}(-1/2)=y^{\prime}(+1/2).
\]

We have $\lambda_{0}=0$ and for $n=0,1,2,\ldots$
\[
\lambda_{4n+1}=(2m\pi+\alpha)^{2};\ \ \lambda_{4n+2}=(2(n+1)\pi-\alpha))^{2};
\]%
\[
\ \lambda_{4n+3}\ =\ \ \lambda_{4n+4}=(2(n+1)\pi))^{2}.
\]
where $\alpha\,=\,\cos^{-1}(-7/8)$

\textbf{Semi-periodic boundary conditions on }$($\textbf{ }$-1/2,+1/2)$%
\textbf{, }\textit{i.e.}%
\[
y(-1/2)=-y(+1/2)\quad\quad y^{\prime}(-1/2)=-y^{\prime}(+1/2).
\]

With $\beta=\cos^{-1}((1+\sqrt{(}33))/16)$ and $\gamma=\cos^{-1}((1-\sqrt
{(}33))/16)$ these are all simple and given by, for $n=0,1,2,\ldots$
\[
\lambda_{4n}\,=\,(2n\pi\,+\,\beta)^{2};\ \lambda_{4n+1}\,=\,(2n\pi
\,+\,\gamma)^{2};\ \
\]%
\[
\lambda_{4n+2}\,=\,(2(n+1)\pi\,-\,\gamma)^{2};\text{ }\lambda_{4n+3}%
\,=\,(2(n+1)\pi\,-\,\beta)^{2}.
\]
See \cite{E} and \cite{Hoc}.

\item \textbf{The Lohner equation}%
\[
-y^{\prime\prime}(x)-1000xy(x)=\lambda y(x)\;\text{for all}\;x\in
(-\infty,+\infty)
\]

Endpoint classification in the space $L^{2}(-\infty,+\infty)$:%
\[%
\begin{tabular}
[c]{cc}\hline
Endpoint & Classification\\\hline
$-\infty$ & LP\\
$+\infty$ & LP\\\hline
\end{tabular}
\]

In \cite{L} Lohner computed the Dirichlet eigenvalues of index (in SLEIGN2
notation) 0, 9, 49 and 99 using interval arithmetic and obtained rigorous
bounds. In double precision SLEIGN2 computed eigenvalues are in good agreement
with these guaranteed bounds.

\item \textbf{The J\"{o}rgens equation}%
\[
-y^{\prime\prime}(x)+(\exp(2x)/4-k\exp(x))y(x)=\lambda y(x)\;\text{for
all}\;x\in(-\infty,+\infty)
\]
where the parameter $k\in(-\infty,+\infty).$

Endpoint classification in the space $L^{2}(-\infty,+\infty),$ for all
$k\in(-\infty,+\infty)$:%
\[%
\begin{tabular}
[c]{cc}\hline
Endpoint & Classification\\\hline
$-\infty$ & LP\\
$+\infty$ & LP\\\hline
\end{tabular}
\]

This is a remarkable example from J\"{o}rgens and SLEIGN2 obtains excellent
results. Details of this problem are given in \cite[Part II, Section 10]{J}.
For all $k\in(-\infty,+\infty)$ the boundary value problem on the interval
$(-\infty,+\infty)$ has a continuous spectrum on $[0,+\infty)$; for $k\leq1/2$
there are no eigenvalues; for $h=0,1,2,3,\ldots$ and then $k$ chosen by
$h<k-1/2\leq h+1,$ there are exactly $h+1$ eigenvalues and these are all below
the continuous spectrum; these eigenvalues are given explicitly by
\[
\lambda_{n}=-(k-1/2-n)^{2},\quad n=0,1,2,3,\ldots,h.
\]

\item \textbf{The Behnke-Goerisch equation}%
\[
-y^{\prime\prime}(x)+k\cos^{2}(x)y(x)=\lambda y(x)\;\text{for all}%
\;x\in(-\infty,+\infty)
\]
where the parameter $k\in(-\infty,+\infty),$

Endpoint classification in the space $L^{2}(-\infty,+\infty),$ for all
$k\in(-\infty,+\infty)$:%
\[%
\begin{tabular}
[c]{cc}\hline
Endpoint & Classification\\\hline
$-\infty$ & LP\\
$+\infty$ & LP\\\hline
\end{tabular}
\]

This is a form of the Mathieu equation. In \cite{BG} these authors computed a
number of Neumann eigenvalues of this problem using interval arithmetic with
rigorous bounds. In double precision SLEIGN2 computed eigenvalues are in good
agreement with these guaranteed bounds.

\item \textbf{The Whittaker equation}%
\[
-y^{\prime\prime}(x)+\left(  \frac{1}{4}+\frac{k^{2}-1}{x^{2}}\right)
y(x)=\lambda\frac{1}{x}y(x)\;\text{for all}\;x\in(0,+\infty)
\]
where the parameter $k\in\lbrack1,+\infty).$

Endpoint classification in the space $L^{2}(0,+\infty),$ for all $k\in
\lbrack1,+\infty)$:%
\[%
\begin{tabular}
[c]{cc}\hline
Endpoint & Classification\\\hline
$0$ & LP\\
$+\infty$ & LP\\\hline
\end{tabular}
\]

This equation is studied in \cite[Part II, Section 10]{J}. There it is shown
that the LP case holds at $+\infty$ and also at $0$ for $k\geq1$. The spectrum
is discrete and is given explicitly by:
\[
\lambda_{n}=n+(k+1)/2,\quad n=0,1,2,3,\ldots.
\]

\item \textbf{The Littlewood-McLeod equation}%
\[
-y^{\prime\prime}(x)+x\sin(x)y(x)=\lambda y(x)\;\text{for all}\;x\in
\lbrack0,+\infty).
\]

Endpoint classification in the space $L^{2}(0,+\infty)$:%
\[%
\begin{tabular}
[c]{cc}\hline
Endpoint & Classification\\\hline
$0$ & R\\
$+\infty$ & LP\\\hline
\end{tabular}
\]

The spectral analysis of this differential equation is considered in
\cite{JEL} and \cite{JBM}; the equation is R at $0$ and LP at $\pm\infty.$ All
self-adjoint operators in $L^{2}[0,\infty)$ have a simple, discrete spectrum
$\{\lambda_{n}:n=0,\pm1,\pm2,\ldots\}$ that is unbounded both above and below,
\textit{i.e.}%
\[
\lim_{n\rightarrow-\infty}\lambda_{n}=-\infty\quad\quad\quad\lim
_{n\rightarrow+\infty}\lambda_{n}=+\infty.
\]
Every eigenfunction has infinitely many zeros in $(0,\infty).$

SLEIGN2, and other codes, fail to compute the eigenvalues for this type of LP
oscillatory problem. However there is qualitative information to be obtained
by considering regular problems on $[0,X]$ with, say, Dirichlet boundary
conditions $y(0)=Y(X)=0.$

\item \textbf{The Morse equation}%
\[
-y^{\prime\prime}(x)+(9\exp(-2x)-18\exp(-x))y(x)=\lambda y(x)\;\text{for
all}\;x\in(-\infty,+\infty)
\]

Endpoint classification in the space $L^{2}(-\infty,+\infty)$%
\[%
\begin{tabular}
[c]{cc}\hline
Endpoint & Classification\\\hline
$-\infty$ & LP\\
$+\infty$ & LP\\\hline
\end{tabular}
\]

This differential equation is studied in \cite[Example 6]{PBB1}; the spectrum
has exactly three negative, simple eigenvalues, and a continuous spectrum on
$[0,\infty);$ the eigenvalues are given explicitly by%
\[
\lambda_{n}=-(n-2.5)^{2}\;\text{for}\;n=0,1,2.
\]

\item \textbf{The Heun equation}%
\[
-(py^{\prime})^{\prime}+qy=\lambda wy\;\text{on}\;(0,1)
\]
where the coefficients $p,q,w$ are given explicitly by, for all $x\in(0,1),$%
\[%
\begin{array}
[c]{l}%
p(x)=x^{c}(1-x)^{d}(x+s)^{e}\\
q(x)=abx^{c}(1-x)^{d-1}(x+s)^{e-1}\\
w(x)=x^{c-1}(1-x)^{d-1}(x+s)^{e-1}.
\end{array}
\]
The parameters $a,b,c,d,e$ and $s$ are all real numbers and satisfy the
following two conditions%
\[
(i)\;s>0\;\text{and}\;c\geq1,d\geq1,a\geq b
\]
and%
\[
(ii)\;a+b+1-c-d-e=0.
\]
From these conditions it follows that%
\[
a\geq1,b\geq1,e\geq1\;\text{and}\;a+b-d\geq1.
\]

The differential equation above is a special case of the general Heun equation%
\[
\frac{d^{2}w(z)}{dz^{2}}+\left(  \frac{\gamma}{z}+\frac{\delta}{z-1}%
+\frac{\varepsilon}{z-a}\right)  \frac{dw(z)}{dz}+\frac{\alpha\beta
z-q}{z(z-1)(z-a)}w(z)=0
\]
with the general parameters $\alpha,\beta,\gamma,\delta,\varepsilon$ replaced
by the real numbers $a,b,c,d,e,$ $a$ replaced by $-s,$ and $q$ replaced by the
spectral parameter $\lambda.$ For general information concerning the Heun
equation see the compendium \cite{AR1}; for the special form of the Heun
equation considered here, and for the connection with confluence of
singularities and applications, see the recent paper \cite{LS1}.

We note that the coefficients of the Sturm-Liouville differential equation
above satisfy the conditions

\begin{itemize}
\item[$(i)$] $q,w\in C[0,1]$ and $w(x)>0$ for all $x\in(0,1)$

\item[$(ii)$] $p^{-1}\in L_{\text{loc}}^{1}(0,1),p(x)>0$ for all $x\in(0,1)$

\item[$(iii)$] $p^{-1}\notin L^{1}(0,1/2]$ and $p^{-1}\notin L^{1}[1/2,1).$
\end{itemize}

Thus both endpoints $0$ and $1$ are singular for the differential equation.
Analysis shows that the endpoint classification for this equation is%
\[%
\begin{tabular}
[c]{ccc}\hline
Endpoint & Parameter & Classification\\\hline
$0$ & $c\in\lbrack1,2)$ & LCNO\\
$0$ & $c\in\lbrack2,+\infty)$ & LP\\
$1$ & $d\in\lbrack1,2)$ & LCNO\\
$1$ & $d\in\lbrack2,+\infty)$ & LP\\\hline
\end{tabular}
\ \
\]

For the endpoint $0$ and for LCNO cases the boundary condition functions $u,v$
are determined by:%
\[%
\begin{tabular}
[c]{ccc}\hline
Parameter & $u$ & $v$\\\hline
$c=1$ & $1$ & $\ln(x)$\\
$1<c<2$ & $1$ & $x^{1-c}$\\\hline
\end{tabular}
\]

For the endpoint $1$ and for LCNO cases the boundary condition functions $u,v$
are determined by:%
\[%
\begin{tabular}
[c]{ccc}\hline
Parameter & $u$ & $v$\\\hline
$d=1$ & $1$ & $\ln(1-x)$\\
$1<d<2$ & $1$ & $(1-x)^{1-d}$\\\hline
\end{tabular}
\]

Further it may be shown that the spectrum of any self-adjoint problem on
$(0,1),$ with the parameters $a,b,c,d,e$ and $s$ satisfying the above
conditions, and considered in the space $L^{2}((0,1);w)$ with either separated
or coupled boundary conditions, is bounded below and discrete. For the
analytic properties, and proofs of the spectral properties of this Heun
differential equation, see the paper \cite{BEHZ}.

In xamples.f, under Example 32, the user can enter a choice for the parameters
$a,b,c,d,e$ and $s,$ subject to the conditions above being satisfied, and
obtain numerical results for the eigenvalues.
\end{enumerate}

\begin{thebibliography}{9}                                                                                                %

\bibitem {AB}M. Abramowitz and I.A. Stegun, \emph{Handbook of mathematical
functions}, Dover Publications, Inc., New York, 1972.

\bibitem {AEZ}F.V. Atkinson, W.N. Everitt and A. Zettl, \emph{Regularization
of a Sturm-Liouville problem with an interior singularity using
quasi-derivatives} Diff. and Int. Equations 1 (1988), 213-222.

\bibitem {PBB1}P.B. Bailey, \emph{SLEIGN}:\emph{ An eigenvalue-eigenfunction
code for Sturm-Liouville problems.} Report Sand77-2044, Sandia National
Laboratory, New Mexico, USA: 1978.

\bibitem {BEHZ}P.B. Bailey, W.N. Everitt, D.B. Hinton and A. Zettl, \emph{Some
spectral properties of the Heun differential equation.} (Submitted for
publication: 2001.)

\bibitem {BEZ}P.B. Bailey, W.N. Everitt and A. Zettl, \emph{Computing
eigenvalues of singular Sturm-Liouville problems}, Results in Mathematics,
\textbf{20 }(1991), 391-423.

\bibitem {BEZ2}P.B. Bailey, W.N. Everitt and A. Zettl, \emph{Regular and
singular Sturm-Liouville problems with coupled boundary conditions}, Proc.
Royal Soc. Edinburgh (A) \textbf{126} (1996), 505-514.

\bibitem {BEZ3}P.B. Bailey, W.N. Everitt and A. Zettl, \emph{The SLEIGN2
Sturm-Liouville code,} ACM Trans. Math. Software. (To appear.)

\bibitem {BEWS}P.B. Bailey, W.N. Everitt, J. Weidmann and A. Zettl,
\emph{Regular approximation of singular Sturm-Liouville problems}. Results in
Mathematics \textbf{23} (1993), 3-22.

\bibitem {BG}H. Behnke and F. Goerisch, \emph{Inclusions for eigenvalues of
self-adjoint problems}, J. Herzberger (Hrsg.): Topics in Validated
Computation, North Holland Elsevier, Amsterdam, (1994).

\bibitem {B}J.P. Boyd, \emph{Sturm-Liouville eigenvalue problems with an
interior pole}, J. Math. Physics, \textbf{22 }(1981), 1575-1590.

\bibitem {DS}N. Dunford and J.T. Schwartz, \emph{Linear Operators, part II},
Interscience Publishers, New York, 1963.

\bibitem {E}M.S.P. Eastham, \emph{The spectral theory of periodic differential
equations}, Scottish Academic Press, Edinburgh and London, 1973.

\bibitem {EGZ}W.N. Everitt, J. Gunson and A. Zettl, \emph{Some comments on
Sturm-Liouville eigenvalue problems with interior singularities}, J. Appl.
Math. Phys. (ZAMP) \textbf{38} (1987), 813-838.

\bibitem {F}G. Fichera, \emph{Numerical and quantitative analysis}, Pitman
Press, London, 1978.

\bibitem {FP1}C.T. Fulton and S. Pruess. `Mathematical software for
Sturm-Liouville problems.' \emph{NSF Final Report for Grants DMS88 and
DMS88-00839; }1993.

\bibitem {H}M.S. Homer, \emph{Boundary value problems for the Laplace tidal
wave equation}, Proc. Roy. Soc. of London (A) \textbf{428} (1990), 157-180.

\bibitem {Hoc}H. Hochstadt, \emph{A special Hill's equation with discontinuous
coefficients}, Amer. Math. Monthly, \textbf{70 }(1963), 18-26.

\bibitem {J}K. J\"{o}rgens, \emph{Spectral theory of second-order ordinary
differential operators}, Lecture Notes: Series no.2, Matematisk Institut,
Aarhus Universitet, 1962/63.

\bibitem {K}A.M. Krall, \emph{Boundary value problems for an eigenvalue
problem with a singular potential}, J. Diff. Equations, \textbf{45} (1982), 128-138.

\bibitem {LS1}W. Lay and S. Yu. Slavyanov, \emph{Heun's equation with nearby
singularities.} Proc. R. Soc. Lond. A \textbf{455} (1999), 4374-4261.

\bibitem {JEL}J.E. Littlewood, \emph{On linear differential equations of the
second order with a strongly oscillating coefficient of }$y$, J. London Math.
Soc. \textbf{41} (1966), 627-638.

\bibitem {L}R.J. Lohner, \emph{Verified solution of eigenvalue problems in
ordinary differential equations}. (Personal communication, 1995.)

\bibitem {M}M. Marletta, \emph{Numerical tests of the SLEIGN software for
Sturm-Liouville problems}, ACM TOMS, \textbf{17 }(1991), 501-503.

\bibitem {JBM}J.B. McLeod, \emph{Some examples of wildly oscillating
potentials}, J. London Math. Soc. \textbf{43} (1968), 647-654.

\bibitem {PMM}P.M. Morse, \emph{Diatomic molecules according to the wave
mechanics}; II: \emph{Vibration levels. }Phys. Rev. \textbf{34} (1929), 57-61.

\bibitem {NZ}H.-D. Niessen and A. Zettl, \emph{Singular Sturm-Liouville
problems; the Friedrichs extension and comparison of eigenvalues}, Proc.
London Math. Soc. \textbf{64} (1992), 545-578.

\bibitem {JDP}J.D. Pryce. \emph{A test package for Sturm-Liouville solvers.}
ACM Trans. Math. Software \textbf{25} (1999), 21-57.

\bibitem {P}M. Plum, \emph{Eigenvalue inclusions for second-order ordinary
differential operators by a numerical homotopy method}, ZAMP, \textbf{41
}(1990), 205-226.

\bibitem {AR1}A. Ronveaux, \emph{Heun differential equations.} (Oxford
University Press; 1995.)

\bibitem {ST}D.B. Sears and E.C. Titchmarsh, \emph{Some eigenfunction
formulae}, Quart. J. Math. Oxford (2) \textbf{1} (1950), 165-175.

\bibitem {T}E.C. Titchmarsh, \emph{Eigenfunction expansions associated with
second- order differential equations}: \textbf{I}. Clarendon Press, Oxford; 1962.

\bibitem {W}G.N. Watson, \emph{A treatise on the theory of Bessel functions},
Cambridge University Press, Cambridge, England, 1958.

\bibitem {Z}A. Zettl, \emph{Computing continuous spectrum}, Proc.
International Symposium, Trends and Developments in Ordinary Differential
Equations, Y. Alavi and P. Hsieh editors, World Scientific, (1994), 393-406.
\end{thebibliography}
\end{document}
SHAR_EOF
fi # end of overwriting check
cd ..
if test ! -d 'Fortran'
then
	mkdir 'Fortran'
fi
cd 'Fortran'
if test ! -d 'Sp'
then
	mkdir 'Sp'
fi
cd 'Sp'
if test ! -d 'Drivers'
then
	mkdir 'Drivers'
fi
cd 'Drivers'
if test -f 'data2'
then
	echo shar: will not over-write existing file "'data2'"
else
cat << "SHAR_EOF" > 'data2'
(output: again: )
(np: param: a: b: classa: classb: bca: bcb: bcc: numeig: end:)
(bca: d n a1,a2)
(bcc:  p s alfa k11,k12,k21,k22)
 
output: results.txt
np: 1 
a: -1.0 
classa:  lcno 
b: 1.0 
classb: lcno 
bca: 1,0  
bcb: 1,0
numeig: 0,2 

again:
a: 0.0
classa: r
bca: d
bcb: 0,1
 
again:

 
np: 2 
param: 0.75
a: 0.0 
classa:  lcno 
b: 1.0 
classb: r 
bca: 1,0  
bcb: d 
numeig: 0,2 

again:
bca: 0,1
 
again:
b: null
bcb: null
classb: lp

again:

np: 3 
param: null
a: 0.0 
classa:  wr 
b: null
classb: lcno 
bca: d  
bcb: 1,0 
numeig: 0,2 

again:

 
np: 4 
a: 0.0 
b: 1.0
classa:  lcno 
classb: r 
bca: 1,0  
bcb: d 
numeig: 0,4 

again:
bca: 0,1
 
again:

 
np: 5 
a: 0.0 
b: 1.0
classa:  wr 
classb: r 
bca: d  
bcb: d 
numeig: 0,4 

again:
bca: 1,1
 
again:
np: 6 
param: null
a: 1.0 
classa:  r 
b: null
classb: lco 
bca: d  
bcb: 1,0 
numeig: -2,2 

again:
bcb: 0,1
 
again:

 
np: 7 
a: 0.0 
b: 1.0
classa: lco 
classb: r 
bca: 1,0  
bcb: d 
numeig: -2,3 

again:
bca: 0,1
numeig: -1,5
 
again:
bca: null
bcb: null
bcc: p
numeig: -1,2

again:
 
bcc: null
param: 1.0
np: 8 
a: 0.0 
b: 1.0
classa:  lcno 
classb: r 
bca: 1,0  
bcb: d 
numeig: 0,1 

again:
bca: 0,1
numeig: 0,2 

again:
param: -0.5
bca: 1,0
numeig: 0,1
 
again:
bca: 0,1
numeig: 0,2

again:

 
param: null
np: 9 
a: 0.0 
b: 1.0
classa: wr
classb: lcno 
bca: d  
bcb: 1,0 
numeig: 0,4 

again:
bcb: 0,1
numeig: 0,4 

again:

 
np: 10 
a: 0.0 
b: 1.0
classa: wr
classb: r 
bca: d  
bcb: d 
numeig: 0,2 

 
again:
np: 11 
a: 0.0 
bca: null
b: pi 
bcb: null
classa: r
classb: r 
bcc: p 
numeig: 0,2 

again:

 
param: 5.0
np: 12 
a: 0.0 
b: pi 
classa: r
classb: r 
bcc: p 
numeig: 0,6 

again:

 
param: -1.0,2.0 
np: 13 
a: 0.0 
classa: lp 
b: null
classb: lp 
bcc: null
numeig: 0,2 

again:

 
param: null
np: 14 
a: 0.0 
classa: r
bca: 5.0,8.0 
classb: lp 
numeig: 0,2 

again:

 
np: 15 
a: null
classa: lp 
bca: null
classb: lp 
numeig: 0,2 

 
again:

np: 16 
param: -0.5,-1.2
a: -1.0
b: 1.0
classa: lp 
classb: wr 
bcb: d
numeig: 0,2 

again:
param: 0.5,-1.2
classb: lcno
bcb: 1,0

again:
param: 1.2,-1.2
classb: lp
bcb: null

again:
param: .5,-.5
classa: wr
classb: lcno
bca: d
bcb: 1,0
 
again:
param: 1.2,-0.5
classb: lp
bcb: null
 
again:
param: 1.2,0.5
classa: lcno
bca: 1,0
 
again:

 
param: null
np: 17 
a: 0.0
classa: lp 
b: null
classb: lp 
bca: null
bcb: null
numeig: 0,3 

again:

 
np: 18 
param: 0.2,0.6
a: -1.0
b: 1.0 
classa: lcno 
classb: lp 
bca: 1,0
numeig: 0,2 

again:

 
np: 19 
param: 0.5,0.2
a: -1.0
b: 1.0 
classa: lco 
classb: lcno 
bca: 1,0
bcb: 1,0
numeig: -1,2 

again:
param: 0.5,0.6
classb: lp
bcb: null
numeig: 0,3

again:

np: 20 
param: 1.0
a: 0.0
classa: lco 
b: null
classb: lp 
bca: 1,0
bcb: null
numeig: -2,1 

 
again:
param: null
np: 21 
a: -pi 
b: pi 
classa: r 
classb: r 
bca: d
bcb: d 
numeig: 0,4
 
again:
bca: null
bcb: null
bcc: p
numeig: 0,4 

again:
bcc: 0.7853982,2.0000000,1.0000000,1.0000000,1.0000000

again:
 
np: 22 
param: 0.0
a: 0.0
classa: lcno 
b: null
classb: lp 
bca: 1,0 
bcb: null
bcc: null
numeig: 0,2
 
again:
bca: 0,1
numeig: 0,2 

again:
param: -1.2
classa: lp
bca: null
 
again:
param: -0.6
classa: wr
bca: d
 
again:
param: 0.5
classa: lcno
bca: 1,0
 
again:
bca: 0,1

again:
param: 1.2
classa: lp
bca: null
 
again:

 
np: 23 
param: 0.0 
a: 0.0
classa: lcno 
b: null
classb: lp 
bca: 1,0 
bcb: null
numeig: 1,2
 
again:
param: -1.2
classa: lp
bca: null
numeig: 0,2 

again:
param: 0.6 
classa: lcno 
bca: 1,0
 
again:
param: 0.5
classa: r
bca: d
 
again:
param: 1.2
classa: lp
bca: null

again:

 
np: 24 
param: -0.6,-1.2
a: -pi/2
b:  pi/2 
classa: lp 
bca: null
classb: lcno 
bcb: 1,0 
numeig: 0,2
 
again:
param: 0.5,-1.2
classb: r 
bcb: d
numeig: 0,2 

again:
param: 1.2,-1.2 
classb: lp 
bcb: null

again:
param: 0.5,-0.6
classa: lcno 
classb: r
bca: 1,0
bcb: d 
 
again:
param: 1.2,-0.6
classb: lp
bca: 1,0
bcb: null

again:
param: 1.2,0.5
classa: r
bca: d

again:

 
param: null
np: 25 
a: -0.5 
b:  0.5 
classa: r 
classb: r 
bca: null
bcb: null
bcc: p
numeig: 0,6

 
again:
param: null
np: 26 
a: 0.0 
b:  1.0 
classa: r 
classb: r 
bca: d
bcb: d
bcc: null
numeig: 0,4

again:
 
np: 27
param: 2.0
a: null
b: null
classa: lp
classb: lp
numeig: 0,2

again:

np: 28 
param: 2.0
a: 0.0
b: pi
classa: r 
classb: r 
bca: n
bcb: n
numeig: 0,6

again:

 
np: 29 
param: 1.0
a: 0.0
classa: lp 
b: null
classb: lp 
bca: null
bcb: null
numeig: 0,2

again:

param: null
np: 30 
a: 0.0
classa: r
b: 10.0 
classb: r 
bca: d 
bcb: d 
numeig: 0,9
 
again:

b: 20.0

again:
b: 30.0

again:
b: 40.0

again:
param: null
np: 31
a: null
b: null
classa: lp
classb: lp
bca: null
bcb: null
numeig: 0,4

again:
np: 32
param: 1.0,4.0,2.0,1.5,4.5
a: 0.0
classa: lcno
bca: 1.0,0.0
b: 1.0
classb: lp
numeig: 0,2

again:
numeig: 18

end:
 


SHAR_EOF
fi # end of overwriting check
if test -f 'driver1.f'
then
	echo shar: will not over-write existing file "'driver1.f'"
else
cat << "SHAR_EOF" > 'driver1.f'
C  PROGRAM COUPDR

C     **********
C     MARCH 1, 2001; P.B. BAILEY, W.N. EVERITT AND A. ZETTL
C     VERSION 1.2
C     **********

      PROGRAM COUPDR

C     This program is for the purpose of indicating how SLCOUP
C     can be called for the purpose of obtaining eigenvalues
C     of Sturm-Liouville problems with regular and singular
C     coupled boundary conditions.
C
C     The call is of the form:
C
C     CALL SLCOUP(A, B, INTAB, P0ATA, QFATA, P0ATB, QFATB,
C    1    A1, A2, B1, B2, NUMEIG, EIG, TOL, IFLAG,
C    2    CPFUN, NCA, NCB, ALFA, K11, K12, K21, K22)
C
C     DECLARE THE NEEDED VARIABLES:

C
C     DECLARE THE NEEDED EXTERNALS:

C
C     SET THE PARAMETERS DEFINING THE INTERVAL OF
C       DEFINITION OF THE DIFFERENTIAL EQUATION:
C       (BESSEL'S EQUATION WITH NU = 0.75)
C
C     .. Scalars in Common ..
      REAL A1S,A2S,AA,ASAV,B1S,B2S,BB,BSAV,DTHDAA,DTHDBB,
     +                 EIGSAV,EPSMIN,FA,FB,GQA,GQB,
     +                 GWA,GWB,HPI,LPQA,
     +                 LPQB,LPWA,LPWB,P0ATAS,P0ATBS,PI,QFATAS,QFATBS,
     +                 TMID,TSAVEL,TSAVER,TWOPI,Z
      INTEGER IND,INTSAV,ISAVE,MDTHZ,MFS,MLS,MMWD,T21
      LOGICAL ADDD,PR
C     ..
C     .. Arrays in Common ..
      REAL TEE(100),TT(7,2),YS(200),YY(7,3,2),ZEE(100)
      INTEGER JAY(100),MMW(100),NT(2)
C     ..
C     .. Local Scalars ..
      REAL A,A1,A2,ALFA,B,B1,B2,EIG,K11,K12,K21,K22,P0ATA,
     +                 P0ATB,QFATA,QFATB,TOL
      INTEGER I,ICPFUN,IFLAG,INTAB,NCA,NCB,NP,NUMEIG
C     ..
C     .. Local Arrays ..
      REAL CPFUN(100),X(100)
C     ..
C     .. External Subroutines ..
      EXTERNAL SLCOUP
C     ..
C     .. Common blocks ..
C     COMMON /ALBE/LPWA,LPQA,LPWB,LPQB
      COMMON /ALBE/LPWA,LPQA,FA,GWA,GQA,LPWB,LPQB,FB,GWB,GQB
      COMMON /BCDATA/A1S,A2S,P0ATAS,QFATAS,B1S,B2S,P0ATBS,QFATBS
      COMMON /DATADT/ASAV,BSAV,INTSAV
      COMMON /DATAF/EIGSAV,IND
      COMMON /LP/MFS,MLS
      COMMON /PASS/YS,MMW,MMWD
      COMMON /PIE/PI,TWOPI,HPI
      COMMON /PRIN/PR,T21
      COMMON /RNDOFF/EPSMIN
      COMMON /TDATA/AA,TMID,BB,DTHDAA,DTHDBB,MDTHZ,ADDD
      COMMON /TEEZ/TEE
      COMMON /TEMP/TT,YY,NT
      COMMON /TSAVE/TSAVEL,TSAVER,ISAVE
      COMMON /Z1/Z
      COMMON /ZEEZ/JAY,ZEE
C     ..
      T21 = 21
      OPEN (T21,FILE='test.out')
      PR = .true.
      A = 0.0D0
      B = 1.0D0

C  ----------------------------------------------------------------C
C     SET THE PARAMETERS DEFINING THE COUPLED BOUNDARY CONDITIONS:
C  ----------------------------------------------------------------C

      A1 = 1.0D0
      A2 = 0.0D0
      B1 = 1.0D0
      B2 = 0.0D0

C     SET THE REMAINING PARAMETERS NEEDED BY SLCOUP:

      P0ATA = -1.0D0
      QFATA = -1.0D0
      P0ATB = -1.0D0
      QFATB = 1.0D0
      INTAB = 1
      NUMEIG = 3
      EIG = 0.0D0
      TOL = 1.D-5
      NCA = 3
      NCB = 1
      ICPFUN = 0

C  ----------------------------------------------------------------C
C     SET THE COUPLED BOUNDARY CONDITION PARAMETERS:
C     (THESE ARE THE "PERIODIC" CONDITIONS.)
C  ----------------------------------------------------------------C
      ALFA = 0.0D0
      K11 = 1.0D0
      K12 = 0.0D0
      K21 = 0.0D0
      K22 = 1.0D0

C ---------------------------------------------------------C
C  If eigenfunction values are wanted, then:

      NP = 42
      DO 10 I = 2,NP - 1
          X(I-1) = A + (I-1)* (B-A)/ (NP-1)
   10 CONTINUE
      ICPFUN = NP - 2

      DO 15 I = 1,ICPFUN
          CPFUN(I) = X(I)
   15 CONTINUE
C  ----------------------------------------------------------------C
C     CALL SLCOUP:
C  ----------------------------------------------------------------C

      CALL SLCOUP(A,B,INTAB,P0ATA,QFATA,P0ATB,QFATB,A1,A2,B1,B2,NUMEIG,
     +            EIG,TOL,IFLAG,ICPFUN,CPFUN,NCA,NCB,ALFA,K11,K12,K21,
     +            K22)

C     WRITE OUT THE RETURNED EIGENVALUE, IT'S ESTIMATED ACCURACY,
C     AND THE IFLAG STATUS:

      WRITE (*,FMT=*) ' NUMEIG, EIG, TOL, IFLAG = ',NUMEIG,EIG,TOL,IFLAG

C
      IF (ICPFUN.GT.0) THEN
          WRITE (*,FMT=*) ' EIGENFUNCTION '
          DO 30 I = 1,ICPFUN
              WRITE (*,FMT=*) X(I),CPFUN(I)
   30     CONTINUE
      END IF
C
      CLOSE (T21)
      STOP
      END
C
      REAL FUNCTION P(X)

C     .. Scalar Arguments ..
      REAL X
C     ..
      P = 1.0D0
      RETURN
      END
C
      REAL FUNCTION Q(X)

C     .. Scalar Arguments ..
      REAL X
C     ..
C     .. Local Scalars ..
      REAL NU
C     ..
      NU = 0.75D0
      Q = (NU*NU-0.25D0)/X**2
      RETURN
      END
C
      REAL FUNCTION W(X)

C     .. Scalar Arguments ..
      REAL X
C     ..
      W = 1.0D0
      RETURN
      END
C
      SUBROUTINE UV(X,U,PUP,V,PVP,HU,HV)

C     .. Scalar Arguments ..
      REAL HU,HV,PUP,PVP,U,V,X
C     ..
C     .. Local Scalars ..
      REAL NU
C     ..
      NU = 0.75D0
      U = X** (NU+0.5D0)
      PUP = (NU+0.5D0)*X** (NU-0.5D0)
      V = X** (-NU+0.5D0)
      PVP = (-NU+0.5D0)*X** (-NU-0.5D0)
      HU = 0.0D0
      HV = 0.0D0
      RETURN
      END
SHAR_EOF
fi # end of overwriting check
if test -f 'driver2.f'
then
	echo shar: will not over-write existing file "'driver2.f'"
else
cat << "SHAR_EOF" > 'driver2.f'
C  PROGRAM DRIVE

C     **********
C     MARCH 1, 2001; P.B. BAILEY, W.N. EVERITT AND A. ZETTL
C     VERSION 1.2
C     **********

      PROGRAM DRIVE
C
C     THIS PROGRAM IS AN "INTERACTIV