\documentclass[reqno]{amsart}
\usepackage{graphicx}
\usepackage{amscd}
\usepackage{subfigure}
\usepackage[arrow,web,tips]{xy}
\usepackage{hyperref}

\AtBeginDocument{{\noindent\small
\emph{Electronic Journal of Differential Equations},
Vol. 2009(2009), No. 15, pp. 1--21.\newline
ISSN: 1072-6691. URL: http://ejde.math.txstate.edu or http://ejde.math.unt.edu
\newline ftp ejde.math.txstate.edu  (login: ftp)}
\thanks{\copyright 2009 Texas State University - San Marcos.}
\vspace{9mm}}

\begin{document}
\title[\hfilneg EJDE-2009/15\hfil Structured matrix numerical solution]
{Structured matrix numerical solution of the nonlinear schr\"odinger
 equation by the inverse scattering transform}

\author[A. Aric\`o, C. van der Mee, S. Seatzu\hfil EJDE-2009/15\hfilneg]
{Antonio Aric\`o, Cornelis van der Mee, Sebastiano Seatzu}

\address{ Dip. Matematica e Informatica,
 Universit\`a di Cagliari, Viale Merello 92, 09123 Cagliari, Italy}
\email[Antonio Aric\`o]{arico@unica.it}
\email[Cornelis van der Mee]{cornelis@krein.unica.it}
\email[Sebastiano Seatzu]{seatzu@unica.it}

\thanks{Submitted December 23, 2008. Published January 13, 2009.}
\thanks{Supported by the Italian  Ministery of Education,
Universities and Research (MIUR) \hfill\break\indent
under  PRIN grant no. 2006017542-003, and by INdAM-GNCS}
\subjclass[2000]{35Q55, 65M32, 45D05}
\keywords{Nonlinear Schr\"odinger  equation; inverse scattering
 transform; \hfill\break\indent
 structured matrices; numerical computation}

\begin{abstract}
  The initial-value problem for the focusing nonlinear
  Schr\"odinger (NLS) equation is solved numerically by
  following the various steps of the inverse scattering
  transform. Starting from the initial value of the solution, a
  Volterra integral equation is solved followed by three FFT to arrive
  at the reflection coefficient and initial Marchenko
  kernel. By convolution these initial data are propagated in
  time. Using structured-matrix techniques the time evolved
  Marchenko integral equation is solved to arrive at the
  solution to the NLS equation.
\end{abstract}

\maketitle
\numberwithin{equation}{section}
\newtheorem{theorem}{Theorem}[section]
\newtheorem{corollary}[theorem]{Corollary}
\newtheorem{lemma}[theorem]{Lemma}
\newtheorem{proposition}[theorem]{Proposition}

\section{Introduction}\label{sec:1}

In the focusing case, the cubic nonlinear Schr\"odinger (NLS)
equation
\begin{equation}\label{1.1}
  \mathbf{i} q_t=q_{xx}+2|q|^2q,\quad q=q(x,t),
\end{equation}
where the subscripts $x\in\mathbb{R}$ and $t>0$ denote partial derivatives,
has many applications \cite{ZS, SCM, AKNS, AS81, NMPZ, AC, APT}. As
usual, we assume the initial potential $q(x,0)$ is
known. This equation arises in signal propagation in optical fibers
\cite{Ag01, HM02, HT73, Sh04}, wave propagation in nonlinear media
\cite{ZS, W74, DEGM}, and evolution of surface waves on sufficiently
deep water \cite{Z68, YL75, LYRF}.

Among the earlier methods to solve the NLS equation numerically are
the split-step and Fourier methods used by Lake et al. \cite{LYRF} and
by Hardin and Tappert \cite{HaT73}. Taha and Ablowitz \cite{TA84a,
TA84b} made a comparative analysis of numerical methods to solve the
NLS equation such as various finite-difference schemes (e.g. the
Ablowitz-Ladik scheme \cite{AL76, AL77}), split-step Fourier methods
\cite{HaT73, LYRF}, and pseudospectral methods \cite{FW73, F96}, and
decided in favor of the split-step Fourier method with pseudospectral
methods as a good second. The dominant numerical method became the
split-step Fourier method (e.g., \cite{WH86, Ag01,XT03}), even though
in the 1990's orthogonal spline collocation methods were successfully
applied to the NLS equation \cite{RFH93, RF94, R97}.

In this paper we propose a method to solve the NLS equation
\eqref{1.1}, based on a new formulation of the inverse scattering
transform (IST) method. As we shall prove, this revised formulation
will allow us to compute the scattering data without resorting to the
Zakharov-Shabat system. The IST method was first developed to solve
the Korteweg-de Vries (KdV) equation by using the direct and inverse
scattering theory of the Schr\"o\-din\-ger equation \cite{GGKM67,
  GGKM74} and subsequently was applied to the NLS equation \cite{ZS}
and other nonlinear evolution equations \cite{AC, AKNS, APT, AS81,
  M86, NMPZ}. Contrary to the numerical methods mentioned above, our
method follows the itinerary of the IST and hence guarantees that the
NLS solutions obtained are indeed the ones found by using the IST method.

 For later use we give some definitions. By $\mathbb{C}^+$
and $\mathbb{C}^-$ we denote the open upper and lower complex half-planes and by
$\overline{\mathbb{C}^+}$ and $\overline{\mathbb{C}^-}$ their closures. Next, let
$H^s(\mathbb{R})$ denote the Sobolev space of those measurable functions $u$ on
$\mathbb{R}$ whose Fourier transforms $\hat{u}$ satisfy
$$
\|u\|_{H^s(\mathbb{R})}:=\Big[\int_{-\infty}^\infty
d\xi\,|\hat{u}(\xi)|^2(1+\xi^2)^s\Big]^{1/2}<+\infty.
$$
Then, apart from a factor $(2\pi)^{-1/2}$, the Fourier transform $\mathcal{F}$ is
a unitary operator from
$L^{2,s}(\mathbb{R}):=L^2(\mathbb{R};(1+\xi^2)^sd\xi)$ onto
$H^s(\mathbb{R})$.

\section{Inverse Scattering Transform Revisited}\label{sec:2}

The inverse scattering transform method to solve the NLS equation
\eqref{1.1} in the so-called focusing case is based on the spectral
analysis of the Zakharov-Shabat system \cite{ZS}
\begin{equation}\label{2.1}
  -\mathbf{i} J\frac{\partial\psi}{\partial x}(\lambda,x)-V(x)\psi(\lambda,x)=\lambda\psi(\lambda,x),
  \quad x\in\mathbb{R},
\end{equation}
where $\lambda$ is a complex spectral parameter,
$$
J=\mathop{\rm diag}(1,-1)=\begin{bmatrix}\,1&0\\ 0&-1\,\end{bmatrix},\quad
V(x)=\begin{bmatrix}0&\mathbf{i}\,q(x,0)\\ \mathbf{i}\,\overline{q(x,0)}&0\end{bmatrix}
=-V(x)^*,
$$
and the complex potential $q(\cdot,0)\in L^1(\mathbb{R})$ is the initial state of the
potential $q(x,t)$. Here overline denotes complex conjugation and the
asterisk the matrix adjoint (conjugate transpose), while the matrix
potential $V(x)$ is skew-selfadjoint. Introducing, as usual in this
context, the right and left Jost matrices $\Psi(\lambda,x)$ and
$\Phi(\lambda,x)$ as those complex $2\times 2$ matrix solutions of
\eqref{2.1} satisfying, for $\lambda\in\mathbb{R}$, the asymptotic conditions
\begin{subequations}\label{2.2}
  \begin{gather}
    \Psi(\lambda,x)=\mathop{\rm diag}(e^{\mathbf{i}\lambda x},e^{-\mathbf{i}\lambda x})+o(1),\quad x\to+\infty,
    \label{2.2a}\\
    \Phi(\lambda,x)=\mathop{\rm diag}(e^{\mathbf{i}\lambda x},e^{-\mathbf{i}\lambda x})+o(1), \quad x\to-\infty,
    \label{2.2b}
  \end{gather}
\end{subequations}
we have, for $\lambda\in\mathbb{R}$, the proportionality relations
\begin{equation}\label{2.3}
  \Psi(\lambda,x)=\Phi(\lambda,x)A_l(\lambda),\quad\Phi(\lambda,x)=\Psi(\lambda,x)A_r(\lambda),
\end{equation}
where the so-called transition matrices $A_l(\lambda)$ and $A_r(\lambda)$ are unitary
$2\times 2$ matrices with unit determinant which are the inverses of each other.
Thus there exist complex numbers $a(\lambda)$ and $b(\lambda)$ such
that\footnote{In the sequel we shall make use of the fact that any
unitary $2\times 2$ matrix with unit determinant has the form
$\left[\begin{smallmatrix}a&b\\
-\overline{b}&\overline{a}\end{smallmatrix}\right]$, where
$|a|^2+|b|^2=1$.}
$$
A_l(\lambda)=\begin{bmatrix}a(\lambda)&b(\lambda)\\ -\overline{b(\lambda)}\;&\;\overline{a(\lambda)}
\end{bmatrix},\quad
A_r(\lambda)=\begin{bmatrix}\overline{a(\lambda)}\;&\;-b(\lambda)\\
\overline{b(\lambda)}&a(\lambda)\end{bmatrix},\quad
 |a(\lambda)|^2+|b(\lambda)|^2=1.
$$
It is easily verified \cite{AKV00} that, for each $\lambda\in\mathbb{R}$, $\Psi(\lambda,x)$
and $\Phi(\lambda,x)$ are unitary matrices with unit determinant.

The direct and inverse scattering theory of the Zakharov-Shabat system
\eqref{2.1} can be found in many sources (e.g., \cite{ZS, SCM, AKNS, EH, AC,
NMPZ, APT}). The direct scattering problem consists of evaluating from
an initial potential $q(x,0)$ in $L^1(\mathbb{R})$ the reflection coefficient
$R(\lambda):=-b(\lambda)/a(\lambda)$ and the
nontrivial so-called bound state solutions of the Zakharov-Shabat
system whose two components belong to $L^2(\mathbb{R})$. The bound state
spectral parameters turn out to be the zeros of the analytic
continuation of $a(\lambda)$ to $\mathbb{C}^+$
and their complex conjugates \cite{AKV00, APT}. On the other hand,
the inverse scattering problem consists of determining the potential
$q(x,t)$, with $q(\cdot,t)\in L^1(\mathbb{R})$ for each $t>0$, from the
reflection coefficient $R(\lambda)$, the bound state spectral parameters,
and finitely many so-called norming constants. These data are used to
write down the coupled system of two Marchenko integral equations
whose unique solution specifies the potential $q(x,t)$ and its energy density
$|q(x,t)|^2$.

Traditionally the inverse scattering transform (IST) to solve the initial-value
problem for the NLS equation \eqref{1.1} is described in terms of the following
three steps:
\begin{itemize}
\item[(1)] Solve the direct scattering problem for \eqref{2.1} with
initial potential $q(x,0)$. Its solution
determines the functions $a(\lambda)$ and $b(\lambda)$ resulting in
the reflection coefficient $R(\lambda)$, the distinct bound state spectral
parameters $\lambda_1,\ldots,\lambda_N\in\mathbb{C}^+$, the norming constants
$\Gamma_{js}$ ($j=1,\ldots,N$, $s=0,1,\ldots,n_j-1$, $n_j$ being the
order of $\lambda_j$ as a zero of $a(\lambda)$), and the initial kernel
$\Omega(\alpha;0)$ of the Marchenko integral equation.

\item[(2)] Propagate the scattering data in time (cf. \cite{ADV07, B})
to arrive at
$$
R(\lambda;t)=e^{4\mathbf{i}\lambda^2t}R(\lambda),\quad\lambda_j(t)=\lambda_j,
$$
$$
\begin{bmatrix}\Gamma_{j0}(t)\\ \Gamma_{j1}(t)\\ \vdots\\ \Gamma_{j,n_j-1}(t)\end{bmatrix}
=e^{-4\mathbf{i} A_j^2t}\begin{bmatrix}\Gamma_{j0}\\ \Gamma_{j1}\\ \vdots\\ \Gamma_{j,n_j-1}
\end{bmatrix},
$$
where the $k,l$-element of the matrix $A_j$ is given by
$\mathbf{i}\lambda_j\delta_{k,l}+\delta_{k,l+1}$ ($k,l=1,\ldots,n_j$), with $\delta_{k,l}$
denoting the Kronecker delta. In particular,
$\Gamma_{j0}(t)=e^{4\mathbf{i}\lambda_j^2t}\Gamma_{j0}$ if $n_j=1$. The propagation of
the Marchenko kernel in time is then computed as a consequence.

\item[(3)] Solve the inverse scattering problem for the time evolved
  scattering data to arrive at the solution $q(x,t)$. This means
  solving the coupled system of Marchenko integral equations,
  associated to the scattering data above.
\end{itemize}

The Marchenko integral equations consist of the following linear system:
\begin{subequations}\label{2.4}
  \begin{gather}
 B_1(x,\alpha;t)=\int_0^\infty d\beta\,\overline{\Omega(2x+\alpha+\beta;t)}\,B_2(x,\beta;t),
    \label{2.4a}\\
 B_2(x,\alpha;t)=-\int_0^\infty d\beta\,\Omega(2x+\alpha+\beta;t)B_1(x,\beta;t)-\Omega(2x+\alpha;t),
    \label{2.4b}
  \end{gather}
\end{subequations}
valid for $\alpha\ge0$, $t\ge0$ and $x\in\mathbb{R}$. The Marchenko kernel $\Omega$ is given by
\begin{equation}\label{2.5}
  \Omega(\alpha;t)=\hat{R}(\alpha;t)+\sum_{j=1}^N\sum_{s=0}^{n_j-1}\,
  \Gamma_{js}(t)\frac{\alpha^s}{s!}e^{\mathbf{i}\lambda_j \alpha}
\end{equation}
for $\alpha\in\mathbb{R}$ and $t\ge0$, where the caret denotes the Fourier transform and
$\hat{R}$ is given
in terms of the reflection coefficient $R(\lambda)$ by the Fourier transform
\begin{equation}\label{2.6}
  \hat{R}(\alpha;t)=\frac{1}{2\pi}\int_{-\infty}^\infty d\lambda\,e^{\mathbf{i}\lambda\alpha}R(\lambda)
  e^{4\mathbf{i}\lambda^2t}.
\end{equation}
As we shall see, the functions $B_1$ and $B_2$, which are the unknown
functions in \eqref{2.4}, are strictly connected to the
Zakharov-Shabat system. The potential $q(x,t)$ and the associated
energy then follow from
\begin{gather}
  q(x,t)=2B_2(x,0^+;t),\label{2.7}\\
  \int_x^\infty dy\,|q(y,t)|^2=-2B_1(x,0^+;t).\label{2.8}
\end{gather}

In this article we do not follow exactly the path laid down by the
inverse scattering transform to solve the initial value problem for
the NLS equation \eqref{1.1}. In fact, our method allows us to obtain the
reflection coefficient $R(\lambda)$ without solving the Zakharov-Shabat system.
%-------------------------------------------------------------------------------
\section{Detailed description of the Method}\label{sec:3}
%-------------------------------------------------------------------------------
Our method for solving the initial value problem for the NLS equation
consists of the following four steps, where the first of three steps
above has been subdivided in two:
\begin{itemize}
  \item[(1)] Compute $B_1(x,\alpha;0)$ and $B_2(x,\alpha;0)$ from the initial data
    $q(x,0)$ by solving the Volterra integral system \eqref{2.4}.

  \item[(2)] The scattering functions $\hat{b}(\lambda)$ and
    $\hat{a}(\lambda)$ are first computed using \eqref{3.14} below
    in terms of $B_1(x,\alpha;0)$ and $B_2(x,\alpha;0)$. Their
    Fourier transforms $b(\lambda)$ and $a(\lambda)$ are then used to evaluate
    the reflection coefficient $R(\lambda)=-b(\lambda)/a(\lambda)$, and determine its
    initial Fourier transform $\hat{R}(\alpha;0)$ from \eqref{2.6}. All
    these steps can be implemented using various FFT. By going through
    the same steps for scattering data involving two different
    truncations of the initial potential $q(x,0)$, we arrive at the
    norming constants $\Gamma_{js}$. In this way
    we solve the direct scattering problem, that is we compute the
    relevant scattering data, without solving the Zakharov-Shabat
    system.

  \item[(3)] The time evolution of $\hat{R}(\alpha;t)$ is then obtained
    by solving the initial-value problem for the linear partial
    differential equation
    \begin{equation}\label{3.1}
\frac{\partial\hat{R}}{\partial t}
+4\mathbf{i}\frac{\partial^2\hat{R}}{\partial\alpha^2}=0
    \end{equation}
    where the initial value $\hat{R}(\alpha;0)$ can be immediately inferred by
    \eqref{2.6}. Adding the norming constant terms as in \eqref{2.5},
    we get the time evolution of the Marchenko kernel $\Omega(\alpha;t)$.

  \item[(4)] Using $\Omega(\alpha;t)$ the Marchenko system \eqref{2.4} is
    then solved, while the solution $q(x,t)$ is evaluated from
    \eqref{2.7} and the associated energy $|q(x,t)|^2$ from
    \eqref{2.8}.
\end{itemize}

Let us organize our method for solving the initial value problem
for the NLS equation in the following six modules.

\subsection{From initial potential $q(x,0)$ to $B_1(x,\alpha)$ and
$B_2(x,\alpha)$}
 From \eqref{2.1} and \eqref{2.2} we easily derive the Volterra
integral equations
\begin{subequations}\label{3.2}
\begin{gather}
\Psi(\lambda,x)=e^{\mathbf{i}\lambda Jx}-\mathbf{i} J\int_x^\infty dy\,e^{-\mathbf{i}\lambda(y-x)J}V(y)
\Psi(\lambda,y),\label{3.2a}\\
\Phi(\lambda,x)=e^{\mathbf{i}\lambda Jx}+\mathbf{i} J\int_{-\infty}^x dy\,e^{-\mathbf{i}\lambda(y-x)J}V(y)
\Phi(\lambda,y),\label{3.2b}
\end{gather}
\end{subequations}
where the entries $\mathbf{i}\,q(x,0)$ and $\mathbf{i}\,\overline{q(x,0)}$ of $V(x)$
belong to $L^1(\mathbb{R})$. Decomposing either of \eqref{3.2} into separate
scalar equations, it is easy to see that the first column of
$\Psi(\lambda,x)e^{-\mathbf{i}\lambda Jx}$ and the second column of
$\Phi(\lambda,x)e^{-\mathbf{i}\lambda Jx}$ are continuous in $\lambda\in\overline{\mathbb{C}^+}$,
are analytic in $\lambda\in\mathbb{C}^+$, and tend to a constant column vector as
$|\lambda|\to\infty$ in $\overline{\mathbb{C}^+}$. Similarly, the second column of
$\Psi(\lambda,x)e^{-\mathbf{i}\lambda Jx}$ and the first column of
$\Phi(\lambda,x)e^{-\mathbf{i}\lambda Jx}$ are continuous in $\lambda\in\overline{\mathbb{C}^-}$,
are analytic in $\lambda\in\mathbb{C}^-$, and tend to a constant column vector as
$|\lambda|\to\infty$ in $\overline{\mathbb{C}^-}$. For details we refer to
\cite{APT, V04, D07}. Now write \cite{V04, VAP05, D07}
\begin{subequations}\label{3.3}
\begin{gather}
\Psi(\lambda,x)=e^{\mathbf{i}\lambda Jx}+\int_x^\infty dy\,K(x,y)e^{\mathbf{i}\lambda Jy},\label{3.3a}\\
\Phi(\lambda,x)=e^{\mathbf{i}\lambda Jx}+\int_{-\infty}^x dy\,G(x,y)e^{\mathbf{i}\lambda Jy},\label{3.3b}
\end{gather}
\end{subequations}
where for each $x\in\mathbb{R}$ the following inequality holds:
$$
\int_x^\infty dy\,|K(x,y)|+\int_{-\infty}^x dy\,|G(x,y)|<+\infty.
$$
We now introduce the matrix functions $B(x,\alpha)$ and $C(x,\alpha)$ by
$$
B(x,\alpha)=K(x,x+\alpha),\quad C(x,\alpha)=G(x,x-\alpha),
$$
where $\alpha\in\mathbb{R}^+$. Then, as functions of $\alpha\in\mathbb{R}^+$, the entries of
$B(x,\alpha)$ and $C(x,\alpha)$ belong to $L^1(\mathbb{R}^+)$. Moreover,
\begin{subequations}\label{3.4}
\begin{gather}
\Psi(\lambda,x)e^{-\mathbf{i}\lambda Jx}=\mathop{\rm diag}(1,-1)+\int_0^\infty d\alpha\,B(x,\alpha)e^{\mathbf{i}\lambda J\alpha},
\label{3.4a}\\
\Phi(\lambda,x)e^{-\mathbf{i}\lambda Jx}=\mathop{\rm diag}(1,-1)+\int_0^\infty d\alpha\,C(x,\alpha)e^{-\mathbf{i}\lambda J\alpha}.
\label{3.4b}
\end{gather}
\end{subequations}
Using the partitioning (see also \cite{AKV00, V04, D07})
$$Z=\begin{bmatrix}Z_1&Z_2\\ Z_3&Z_4\end{bmatrix},$$
we obtain, from the unitarity of $\Psi(\lambda,x)$ and $\Phi(\lambda,x)$, that
\begin{subequations}\label{3.5}
\begin{align}
B_1(x,\alpha)&=\overline{B_4(x,\alpha)},\quad B_2(x,\alpha)=-\overline{B_3(x,\alpha)},
\label{3.5a}\\
C_1(x,\alpha)&=\overline{C_4(x,\alpha)},\quad C_2(x,\alpha)=-\overline{C_3(x,\alpha)}.
\label{3.5b}
\end{align}
\end{subequations}

The following result has been proved in detail in \cite{D07, VAP05}, while in
\cite{AKV00} the proof is given for potentials $V(x)$ satisfying $V(x)^*=V(x)$.
Observe that Theorem \ref{th:3.1} and the second of \eqref{3.5a} imply
\eqref{2.7} and \eqref{2.8} for $t=0$.

\begin{theorem}\label{th:3.1}
We have the following identities:
\begin{subequations}\label{3.6}
\begin{gather}
B_1(x,\alpha)=\int_x^\infty dy\,q(y,0)B_3(y,\alpha),\label{3.6a}\\
B_3(x,\alpha)=-\frac{1}{2}\,\overline{q(x+\alpha/2,0)}
-\int_x^{x+\frac{\alpha}{2}}dy\,\overline{q(y,0)}B_1(y,\alpha-2(y-x)).\label{3.6b}
\end{gather}
\end{subequations}
Moreover, \eqref{3.6} has a unique solution satisfying
\begin{equation}\label{3.7}
\sup_{x\in\mathbb{R}}\,\int_0^\infty dy\left(|B_1(x,\alpha)|+|B_3(x,\alpha)|\right)
<+\infty.
\end{equation}
\end{theorem}

Equations \eqref{3.6} suggest a method for computing $B_1(x,\alpha)$ and
$B_2(x,\alpha)$ from the potential $q(x,0)$ when also using the second symmetry
relation \eqref{3.5a}. Alternatively, the linear system
\begin{subequations}\label{3.8}
\begin{gather}
C_1(x,\alpha)=-\int_{-\infty}^x dy\,q(y,0)C_3(y,\alpha),\label{3.8a}\\
C_3(x,\alpha)=\frac{1}{2}\,\overline{q(x-\alpha/2,0)}
+\int_{x-\alpha/2}^x dy\,\overline{u(y)}\,C_1(y,\alpha+2(y-x)),\label{3.8b}
\end{gather}
\end{subequations}
is uniquely solvable under the condition that
$$\sup_{x\in\mathbb{R}}\,\int_0^\infty dy\left(|C_1(x,\alpha)|+|C_3(x,\alpha)|\right)
<+\infty.$$
Using the second of the symmetry relations \eqref{3.5b} we can compute
$C_1(x,\alpha)$ and $C_2(x,\alpha)$ from $q(x,0)$ by solving \eqref{3.8}.

\subsection{From $B_1(x,\alpha)$ and $B_2(x,\alpha)$ to the Marchenko kernel}

Writing the reflection coefficient in the form
\begin{equation}\label{3.9}
R(\lambda)=\int_{-\infty}^\infty d\alpha\,e^{-\mathbf{i}\lambda\alpha}\hat{R}(\alpha),
\end{equation}
the Marchenko integral equations have the form \eqref{2.4}, where, putting
$\Omega(\alpha)=\Omega(\alpha;0)$ and $\hat{R}(\alpha)=\hat{R}(\alpha,0)$, we have
\begin{equation}\label{3.10}
\Omega(\alpha)=\hat{R}(\alpha)+\sum_{j=1}^N\sum_{s=0}^{n_j-1}\,
\Gamma_{js}\frac{\alpha^s}{s!}e^{\mathbf{i}\lambda_j\alpha}.
\end{equation}
Here $\lambda_1,\ldots,\lambda_N$ are the distinct bound state spectral parameters in
$\mathbb{C}^+$ and $\Gamma_{js}$ are the corresponding norming constants, where
$\Gamma_{j,n_j-1}\neq0$ \cite{ADV07, D07}. Recalling that $a(\lambda)$ extends to
a function that is continuous in $\lambda\in\overline{\mathbb{C}^+}$, is analytic in
$\lambda\in\mathbb{C}^+$, and tends to $1$ as $|\lambda|\to\infty$ in $\overline{\mathbb{C}^+}$, we
\textbf{assume}\footnote{A technical assumption of this type appears in
\cite{AKNS, APT, D07}.} that $a(\lambda)\neq0$ for $\lambda\in\mathbb{R}$ and define the
transmission coefficient $T(\lambda)$ by
$$
T(\lambda)=\frac{1}{a(\lambda)}.
$$
Then it easily follows that the bound
state spectral parameter values $\lambda_j\in\mathbb{C}^+$ are the (finitely many)
poles of $T(\lambda)$.

The following theorem shows that for sufficiently large $x\in\mathbb{R}$ the Marchenko
kernel $\Omega(\gamma)=\Omega(\gamma;0)$ can be found from $B_1(x,\alpha)$ and $B_2(x,\alpha)$ by
solving \eqref{2.4b} for $\Omega$ \cite[App.~B.1]{D07}.

\begin{theorem}\label{th:3.2}
Suppose that, for some $x\in\mathbb{R}$, the function
$$
\Psi_1(\lambda,x)=e^{\mathbf{i}\lambda x}\Big[1+\int_0^\infty d\alpha\,e^{\mathbf{i}\lambda\alpha}B_1(x,\alpha)
\Big]=e^{\mathbf{i}\lambda x}+\int_x^\infty dy\,e^{\mathbf{i}\lambda y}K_1(x,y)
$$
does not vanish for $\lambda\in\overline{\mathbb{C}^+}$. Then the convolution integral
equation
\begin{equation}\label{3.11}
\Omega(x;\alpha)+\int_\alpha^\infty d\gamma\,B_1(x,\gamma-\alpha)\Omega(x;\gamma)=-B_2(x,\alpha),\quad
\alpha\in\mathbb{R}^+,
\end{equation}
has a unique solution satisfying
$$\int_0^\infty d\alpha\,|\Omega(x;\alpha)|<+\infty.$$
Moreover, there exists $x_0\in\mathbb{R}$ such that \eqref{3.11} is uniquely
solvable for $x>x_0$.
\end{theorem}

\begin{proof}
If $\Psi_1(\lambda,x)\neq0$ for each $\lambda\in\overline{\mathbb{C}^+}$, there exists
$Z(x,\alpha)$ with $Z(x,\cdot)\in L^1(\mathbb{R}^+)$ such that
$$
\frac{e^{\mathbf{i}\lambda x}}{\Psi_1(\lambda,x)}=1+\int_0^\infty d\alpha\,e^{\mathbf{i}\lambda\alpha}Z(x,\alpha),
\quad\lambda\in\overline{\mathbb{C}^+}.
$$
Using that
$$
Z(x,\alpha)+B_1(x,\alpha)+\int_0^\alpha d\gamma\,B_1(x,\gamma)Z(\alpha-\gamma)=0,\quad\alpha\in\mathbb{R}^+,
$$
it is easily verified that
$$
\Omega(x;\alpha)=-B_2(x,\alpha)-\int_\alpha^\infty d\gamma\,Z(x,\gamma-\alpha)B_2(x,\gamma)
$$
is the unique solution of \eqref{3.11}.
Equations \eqref{3.6a} and \eqref{3.7} in combination with
$q(\cdot,0)\in L^1(\mathbb{R})$
imply that
$$
\lim_{x\to+\infty}\,\int_0^\infty d\alpha\,|B_1(x,\alpha)|=0.
$$
Choosing $x_0\in\mathbb{R}$ such that $\int_0^\infty d\alpha\,|B_1(x,\alpha)|<1$ for
$x>x_0$, we get
$$
|\Psi_1(\lambda,x)-e^{\mathbf{i}\lambda x}|<1,\quad x>x_0\ \text{and}\
\lambda\in\overline{\mathbb{C}^+},
$$
which implies the unique solvability of \eqref{3.11} for $x>x_0$.
\end{proof}

Equations \eqref{2.2b} and \eqref{2.3} imply
$$
\lim_{x\to-\infty}\sup_{\lambda\in\overline{\mathbb{C}^+}}\,
|\Psi_1(\lambda,x)-e^{\mathbf{i}\lambda x}a(\lambda)|=0.
$$
Thus there exists $x_1\in\mathbb{R}$
such that $\Psi_1(\lambda,x)$ does not vanish for $x<x_1$ and
$\lambda\in\overline{\mathbb{C}^+}$, i.e., \eqref{3.11} is uniquely solvable for
$x<x_1$ if and only if $a(\lambda)\neq0$ for every
$\lambda\in\overline{\mathbb{C}^+}$, i.e., if and only if there are no bound
states. By the argument principle \cite{C73} it then easily follows
that the combined zero order of $\Psi_1(\lambda,x)$, under the hypothesis
that $\Psi_1(\lambda,x)\neq0$ for $\lambda\in\mathbb{R}$, is given by
$$
p(x)=\frac{1}{2\pi}\left[\arg\left(e^{-\mathbf{i}\lambda x}\Psi_1(\lambda,x)\right)
\right]_{\lambda=-\infty}^{+\infty}.
$$
We have $p(x)=0$ for $x>x_0$ and $p(x)$ equal to the combined
pole order of the transmission coefficient for $x$ in some left half-line,
with jumps occurring at those values of $x$ where $\Psi_1(\lambda,x)$ has a zero
$\lambda\in\mathbb{R}$.

\begin{corollary}\label{th:3.3}
Suppose $x\in\mathbb{R}$. Then the following statements are true:
\begin{itemize}

\item[(1)] Equation \eqref{3.11} has at least one solution $\Omega(x;\alpha)$
such that
$$\int_0^\infty d\alpha\,|\Omega(x;\alpha)|<+\infty$$
and
\begin{equation}\label{3.12}
\int_\alpha^\infty d\gamma\,\overline{B_2(x,\gamma-\alpha)}\Omega(x;\gamma)
=\overline{B_1(x,\alpha)},\quad\alpha\in\mathbb{R}^+.
\end{equation}

\item[(2)] Equation \eqref{3.11} is a Fredholm integral equation of the
second kind if and only if $\Psi_1(\lambda,x)\neq0$ for $\lambda\in\mathbb{R}$.

\item[(3)] The number of linearly independent solutions of the homogeneous
counterpart of \eqref{3.11} coincides with $p(x)$.
\end{itemize}
\end{corollary}

\begin{proof}
The first part follows by observing that the solution $\Omega(x;\alpha)=\Omega(2x+\alpha)$
of \eqref{3.11} also satisfies \eqref{3.12}, because \eqref{3.11} and \eqref{3.12}
together coincide with the coupled Marchenko system \eqref{2.4} on switching known
and unknown functions. Next, applying the Fourier transform to \eqref{3.11} for
$\alpha\in\mathbb{R}^+$, we obtain
$$e^{\mathbf{i}\lambda x}\Psi_1(-\lambda,x)\hat{\Omega}(x;\lambda)=-\hat{B}_2(x,\lambda),\quad\lambda\in\mathbb{R},$$
which implies the second and third parts.
\end{proof}

\subsection{Separation of the Marchenko kernel}

To implement the time evolution of
the initial Marchenko kernel $\Omega(\alpha)$ without knowing in advance the reflection and
bound state contributions, we proceed as follows. Singling out the first row of
the limit of $\Psi(\lambda,x)e^{-\mathbf{i}\lambda Jx}$ as $x\to-\infty$, we obtain from
\eqref{3.2a} with the help of \eqref{2.2a} and \eqref{2.3}
\begin{gather*}
e^{-\mathbf{i}\lambda x}\Psi_1(\lambda,x)
=1+\int_x^\infty dy\,q(y,0)\Psi_3(\lambda,y)e^{-\mathbf{i}\lambda y},\\
e^{-\mathbf{i}\lambda x}\Psi_2(\lambda,x)=\int_x^\infty dy\,e^{-\mathbf{i}\lambda y}q(y,0)\Psi_4(\lambda,y).
\end{gather*}
Taking the limit as $x\to-\infty$ and using \eqref{3.4a} we get
\begin{subequations}\label{3.13}
\begin{gather}
a(\lambda)=1+\int_0^\infty d\alpha\,e^{\mathbf{i}\lambda\alpha}\hat{a}(\alpha),\label{3.13a}\\
b(\lambda)=\int_{-\infty}^\infty d\alpha\,e^{-\mathbf{i}\lambda\alpha}\hat{b}(\alpha),\label{3.13b}
\end{gather}
\end{subequations}
where [cf. \eqref{3.5a}]
\begin{subequations}\label{3.14}
\begin{equation}
\hat{a}(\alpha)=\int_{-\infty}^\infty dy\,q(y,0)B_3(y,\alpha)
=-\int_{-\infty}^\infty dy\,q(y,0)\,\overline{B_2(y,\alpha)},\label{3.14a}
\end{equation}
\begin{align}
\hat{b}(\alpha)&=\frac{1}{2}\Big[q(\alpha/2,0)+\int_{-\infty}^\alpha dy\,
q(y/2,0)B_4(y/2,\alpha-y)\Big]\nonumber\\
&=\frac{1}{2}\Big[q(\alpha/2,0)+\int_{-\infty}^\alpha dy\,q(y/2,0)\,
\overline{B_1(y/2,\alpha-y)}\Big].\label{3.14b}
\end{align}
\end{subequations}
Using truncated potentials as specified below, the numerical solution of
\eqref{3.14} can be greatly simplified.

\subsection{Computation of the Initial Marchenko Kernel}

 To compute the initial Marchenko kernel without solving the
Za\-kha\-rov-Sha\-bat system, we introduce a technique based on the
use of truncated potentials. For $\xi\in\mathbb{R}$ we define as a truncated
potential the function
\begin{equation}\label{3.15}
q_\xi(x,0)=\begin{cases}q(x,0),&x>\xi,\\ 0,&x<\xi,\end{cases}
\end{equation}
which is associated to the initial potential $q(x,0)$.
Then the corresponding right Jost matrix is given by
$$
\Psi_\xi(\lambda,x)=\begin{cases}\Psi(\lambda,x),&x\ge\xi,\\ e^{\mathbf{i}\lambda J(x-\xi)}
\Psi(\lambda,\xi),&x\le\xi.\end{cases}
$$
Hence
\begin{subequations}\label{3.16}
\begin{gather}
a_\xi(\lambda)=e^{-\mathbf{i}\lambda\xi}\Psi_1(\lambda,\xi),\quad
T_\xi(\lambda)=\frac{1}{a_\xi(\lambda)}=\frac{e^{\mathbf{i}\lambda\xi}}{\Psi_1(\lambda,\xi)},
\label{3.16a}\\
b_\xi(\lambda)=e^{-\mathbf{i}\lambda\xi}\Psi_2(\lambda,\xi),\quad
R_\xi(\lambda)=-\frac{b_\xi(\lambda)}{a_\xi(\lambda)}
=-\frac{\Psi_2(\lambda,\xi)}{\Psi_1(\lambda,\xi)}.\label{3.16b}
\end{gather}
\end{subequations}
Consequently, $p(\xi)$ coincides with the combined pole order of the
transmission coefficient $T_\xi(\lambda)$ corresponding to the truncated
potential $q_\xi$.

Let us now explain an easy algorithm for computing separately the
contributions to the initial Marchenko kernel of the reflection
coefficient and of the bound states. First we compute $B_1(x,\alpha)$ and
$B_2(x,\alpha)$ from the potential $q(x,0)$. Then $\hat{a}(\alpha)$ and
$\hat{b}(\alpha)$ are computed with the help of \eqref{3.14}. These
Fourier transforms are used to get $b(\lambda)$ and $a(\lambda)$, to evaluate
the initial reflection coefficient $R(\lambda)=-b(\lambda)/a(\lambda)$, and
finally to compute its Fourier transform $\hat{R}(\alpha)$.

 We note that an appropriate use of the truncated potential may
 greatly simplify the computation of $\hat{a}(\alpha)$ and
 $\hat{b}(\alpha)$. Indeed, denoting by $\hat{a}_\xi$ and $\hat{b}_\xi$
 the functions $\hat{a}$ and $\hat{b}$ associated with the truncated
 potential, we have
\begin{equation}\label{3.17}
\hat{a}_\xi(\alpha)=\begin{cases}B_1(\xi,\alpha),&\alpha\ge0,\\ 0,&\alpha<0,\end{cases}
\quad
\hat{b}_\xi(\alpha)=\begin{cases}-B_3^*(\xi,\alpha-2\xi),&\alpha\ge2\xi,\\
0,&\alpha<2\xi.\end{cases}
\end{equation}

If $\|q\|_1\le(\pi/2)$, there are no bound states \cite{KS03} and hence
$\Omega(\alpha)=\hat{R}(\alpha)$. On the other hand, if $\|q\|_1>(\pi/2)$, we
choose $\zeta\in\mathbb{R}$ such that $\|q_\zeta\|_1<(\pi/2)$ and compute the
corresponding $\hat{R}_\zeta(\alpha)$. Further, we compute the Marchenko
kernel $\Omega_\zeta(\alpha)$ by solving \eqref{2.4} at $t=0$. Then
the bound state part of the Marchenko kernel coincides with
$\Omega_\zeta(\alpha)-\hat{R}_\zeta(\alpha)$ and it can be identified by a least
squares technique.

\subsection{Time evolution of the Marchenko kernel}

The Fourier transformed
reflection coefficient $\hat{R}(\alpha;t)$ and the Marchenko kernel
$\Omega(\alpha;t)$ both satisfy the partial differential equation \eqref{3.1}.
The initial value problem for this PDE has a unique solution in each Sobolev
space $H^s(\mathbb{R})$ and its solution is represented by a strongly continuous group
of unitary operators $U(t)$ acting on the initial data. This is easily verified
by moving \eqref{3.1} to $L^{2,s}(\mathbb{R})$ by using $(2\pi)^{-1/2}\mathcal{F}$ as a unitary
equivalence, as indicated by the  diagram
$$
\begin{CD}
L^{2,s}(\mathbb{R})@>{e^{4\mathbf{i}\lambda^2t}}>>L^{2,s}(\mathbb{R})\\
@V{(2\pi)^{-1/2}\mathcal{F}}VV @VV{(2\pi)^{-1/2}\mathcal{F}}V\\
H^s(\mathbb{R})@>>{U(t)}>H^s(\mathbb{R})
\end{CD}
$$

In \cite{ADV07} the time evolution of the Marchenko kernel has been
given in the multisoliton case, where the reflection coefficient
$R(\lambda)\equiv0$. In this case there exist a complex $p\times p$ matrix $A$
having all of its eigenvalues in the right half-plane, a complex column vector
$B$ and a complex row vector $C$ such that
$$
\Omega(\alpha;0)=Ce^{-\alpha A}B.
$$
It then follows that the time evolved Marchenko kernel is given by
$$
\Omega(\alpha;t)=Ce^{-\alpha A}e^{-4\mathbf{i} tA^2}B.
$$

The free Schr\"o\-din\-ger group closely related to $U(t)$ has been studied in
detail in \cite[Ch.~2]{C03}, where Lemma \ref{th:3.4} is proved in a different
way.

\begin{lemma}\label{th:3.4}
Let $\hat{R}(\,\cdot\,;0)$ be an element of $L^1(\mathbb{R})\cap L^2(\mathbb{R})$. Then the
unique solution $\hat{R}(\,\cdot\,;t)$ of the partial differential equation
\begin{equation}\label{3.18}
\frac{\partial\hat{R}}{\partial t}
+4\mathbf{i}\frac{\partial^2\hat{R}}{\partial\alpha^2}=0
\end{equation}
in $L^2(\mathbb{R})$ is given by the expression
\begin{equation}\label{3.19}
\hat{R}(\alpha;t)=[U(t)\hat{R}(\,\cdot\,;0)](\alpha)=\frac{1+\mathbf{i}}{4\sqrt{2\pi|t|}}
\int_{-\infty}^\infty d\beta\,e^{-\frac{1}{16t}\mathbf{i}(\alpha-\beta)^2}\hat{R}(\beta;0).
\end{equation}
If $\hat{R}(\,\cdot\,;0)\in L^2(\mathbb{R})$, the solution of \eqref{3.18} in $L^2(\mathbb{R})$
is given by
\begin{equation}\label{3.20}
\hat{R}(\alpha;t)=[U(t)\hat{R}(\,\cdot\,;0)](\alpha)=\lim_{N\to+\infty}\,
\frac{1+\mathbf{i}}{4\sqrt{2\pi|t|}}\int_{-N}^N d\beta\,e^{-\frac{1}{16t}\mathbf{i}(\alpha-\beta)^2}
\hat{R}(\beta;0).
\end{equation}
\end{lemma}

\begin{proof}
Suppose $\hat{R}(\,\cdot\,;0)\!\in\!L^2(\mathbb{R})$. Then using Plancherel's theorem
twice we get
\begin{align*}
\hat{R}(\alpha;t)&=\lim_{N\to+\infty}\frac{1}{2\pi}\int_{-N}^N d\lambda\,
e^{-\mathbf{i}\lambda\alpha} \hat{\hat{R}}(\lambda;t)\!=\!\lim_{N\to+\infty}\frac{1}{2\pi}\int_{-N}^N d\lambda\,
e^{-\mathbf{i}\lambda\alpha}e^{4\mathbf{i}\lambda^2t}\hat{\hat{R}}(\lambda;0)\\ &=\lim_{N\to+\infty}
\frac{1}{2\pi}\int_{-N}^N d\lambda\,e^{-\mathbf{i}\lambda\alpha}e^{4\mathbf{i}\lambda^2t}
\int_{-\infty}^\infty d\beta\,e^{\mathbf{i}\lambda\beta}\hat{R}(\beta;0).
\end{align*}
Supposing that $\hat{R}(\,\cdot\,;0)\!\in\!L^1(\mathbb{R})\cap L^2(\mathbb{R})$, we apply
 Fubini's Theorem to obtain
\begin{align*}
\hat{R}(\alpha;t)
&=\lim_{N\to+\infty}\int_{-\infty}^\infty d\beta\Big(\frac{1}{2\pi}
\int_{-N}^N d\lambda\,e^{-\mathbf{i}\lambda(\alpha-\beta)}e^{4\mathbf{i}\lambda^2t}\Big)\hat{R}(\beta;0)\\
&=\lim_{N\to+\infty}\int_{-\infty}^\infty d\beta\Big(\frac{1}{2\pi}
\int_{-N}^N d\lambda\,e^{4\mathbf{i} t[\lambda-\frac{1}{8t}(\alpha-\beta)]^2}\Big)
e^{-\frac{1}{16t}\mathbf{i}(\alpha-\beta)^2}\hat{R}(\beta;0).
\end{align*}
Making the change of variable $\eta=2\sqrt{|t|}[\lambda-\frac{1}{8t}(\alpha-\beta)]$,
we now get
\begin{equation}\label{3.21}
\hat{R}(\alpha;t)=\lim_{N\to+\infty}\int_{-\infty}^\infty d\beta\Big(
\frac{1}{4\pi\sqrt{|t|}}\int_{2\sqrt{|t|}[-N-\frac{1}{8t}(\alpha-\beta)]}
^{2\sqrt{|t|}[N-\frac{1}{8t}(\alpha-\beta)]}d\eta\,e^{\mathbf{i}\eta^2}\Big).
\end{equation}
Using the Fresnel integrals \cite[7.3.1-7.3.2]{AS64}
$$
C(z)=\int_0^z dt\,\cos\left(\frac{\pi}{2}t^2\right),\quad
S(z)=\int_0^z dt\,\sin\left(\frac{\pi}{2}t^2\right),
$$
which satisfy $C(+\infty)=S(+\infty)=\frac{1}{2}$ (cf. \cite[7.3.0]{AS64}),
we easily calculate
\begin{align*}
&\lim_{N\to+\infty}
\frac{1}{4\pi\sqrt{|t|}}\int_{2\sqrt{|t|}[-N-\frac{1}{8t}(\alpha-\beta)]}
^{2\sqrt{|t|}[N-\frac{1}{8t}(\alpha-\beta)]}d\eta\,e^{\mathbf{i}\eta^2}
=\frac{1}{2\pi\sqrt{|t|}}\int_0^\infty d\eta\,e^{\mathbf{i}\eta^2}\\
&=\frac{1}{2\pi\sqrt{|t|}}\Big[\int_0^\infty d\eta\,\cos(\eta^2)
+\mathbf{i}\int_0^\infty d\eta\,\sin(\eta^2)\Big]=\frac{1+\mathbf{i}}{4\sqrt{2\pi|t|}}.
\end{align*}
Using the Theorem of Dominated Convergence, \eqref{3.21} can be written as
\begin{align*}
\hat{R}(\alpha;t)&=\int_{-\infty}^\infty d\beta\lim_{N\to+\infty}\Big(
\frac{1}{4\pi\sqrt{|t|}}\int_{2\sqrt{|t|}[-N-\frac{1}{8t}(\alpha-\beta)]}
^{2\sqrt{|t|}[N-\frac{1}{8t}(\alpha-\beta)]}d\eta\,e^{\mathbf{i}\eta^2}\Big)\\
&\quad \times e^{-\frac{1}{16t}\mathbf{i}(\alpha-\beta)^2}\hat{R}(\beta;0)\\
&=\frac{1+\mathbf{i}}{4\sqrt{2\pi|t|}}\int_{-\infty}^\infty d\beta\,
e^{-\frac{1}{16t}\mathbf{i}(\alpha-\beta)^2}\hat{R}(\beta;0),
\end{align*}
which implies \eqref{3.19}. Equation \eqref{3.20} now easily follows.
\end{proof}

We now derive the following result.

\begin{theorem}\label{th:3.5}
Suppose $\hat{R}\in L^1(\mathbb{R})\cap L^2(\mathbb{R})$. Then the Marchenko kernel is
given by
\begin{equation}\label{3.22}
\Omega(\alpha;t)=\hat{R}(\alpha;t)+Ce^{-\alpha A}e^{-4\mathbf{i} tA^2}B,
\end{equation}
where $A$, $B$ and $C$ are matrices and
\begin{equation}\label{3.23}
\hat{R}(\alpha;t)=\frac{1+\mathbf{i}}{4\sqrt{2\pi|t|}}\int_{-\infty}^\infty d\beta\,
e^{-\frac{1}{16t}\mathbf{i}(\alpha-\beta)^2}\hat{R}(\beta;0).
\end{equation}
\end{theorem}

\subsection{Solving Marchenko equations to get NLS solution}

Now that the time evolved Marchenko kernel $\Omega(\alpha,t)$ has been evaluated,
we solve the system of Marchenko equations \eqref{2.4} and then apply
\eqref{2.7} to arrive at the NLS solution $q(x,t)$. To compute the
squared modulus $|q(x,t)|^2$, we employ \eqref{2.8} in the following
form:
\begin{equation}\label{3.24}
|q(x,t)|^2=2\frac{\partial B_1}{\partial x}(x,0^+;t).
\end{equation}

\section{Algorithms}\label{sec:4}

To solve the NLS equation \eqref{1.1} we start with a truncated
potential $q_\xi(x,0)$ and apply six steps, under the hypothesis that
$q(x,0)$ has at most one bound state. We note that our method recognizes
automatically if the initial potential allows one bound state or none.
In other words, assume that $q_\xi$ has at most one bound state.

The six steps are as follows:

\subsection{Solution of the Volterra system}

We have to solve the system
\begin{subequations}\label{4.1}
\begin{gather}
B_{1,\xi}(x,\alpha)-\int_x^{+\infty} q_\xi(y) B_{3,\xi}(y,\alpha)d{y}=0,
\label{4.1a}\\
\int_x^{x+\frac{\alpha}{2}}\overline{q_\xi(y)}B_{1,\xi}
\big(y,\alpha-2(y-x)\big)d{y} +B_{3,\xi}(x,\alpha)
=-\frac12\overline{q_\xi(x+\alpha/2)}.\label{4.1b}
\end{gather}
\end{subequations}
where $q_\xi(x,0)$, $x\in\mathbb{R}$, is given and $B_{1,\xi}(x,\alpha)$
and $B_{3,\xi}(x,\alpha)$, $x\in\mathbb{R}$ and $\alpha\in\mathbb{R}^+$, are
the functions to compute. The system is solved in a Cartesian
grid by recursively applying a local formula obtained by the composite
Simpson's rule. The recursion scheme is applied in each row
$\{\alpha_j=2hj\}$ of the grid starting with $j=2$, and going on from
right to left in each row grid $\{x_i=ih\}$, as Fig.~\ref{fig:vols}
shows.
\begin{figure}[ht]
\centering
  \mbox{
  \def\latticebody{\drop{.}}
  \begin{xy}0;<1cm,0cm>:
    %
    \ar@{->}(-1.3,0);(4.3,0)         \POS (4.1,-.2)*{x}
    \ar@{->}(0,-.5);(0,3.6)         \POS (-.2,3.2)*{\alpha}
    %
%    \POS (0,0)*=0{\circ}
    \POS 0;<5mm,0cm>:<0cm,10mm>::
    \xylattice{-1}{7}{-1}{3}%
    %
    \POS 0;<1cm,0cm>:
    \POS (1  ,3)*=0{\bullet}
    \POS (1.5,3)*=0{\circ}
    \POS (2  ,3)*=0{\circ}
    \POS (1.5,2)*=0{\circ}
    \POS (2  ,1)*=0{\circ}
    \ar@{-}(1,3);(2,3)
    \ar@{-}(1,3);(2,1)
    %
    \SelectTips{cm}{10}
    \ar@{..}(3,1);(3,0.8)
    \ar@{..}(3.5,1);(3.5,.8)
    \ar@{<->}@<-2mm>(3,1);(3.5,1)   _{h}
    \ar@{..}(3.5,1);(3.7,1)
    \ar@{..}(3.5,2);(3.7,2)
    \ar@{<->}@<-2mm>(3.5,1);(3.5,2) _{2h}
  \end{xy}}
\caption{Pattern used for solving the Volterra system \eqref{4.1}}
\label{fig:vols}
\end{figure}
Note that $B_{1,\xi}$, $B_{3,\xi}$ can be obtained by using the
relations
\begin{subequations}\label{4.2}
\begin{align}
B_{1,\xi}(x,\alpha)&=\begin{cases}B_1(x,\alpha),&x\geq\xi,\\
B_1(\xi,\alpha),&x<\xi,\end{cases}\label{4.2a}\\
B_{3,\xi}(x,\alpha)&=\begin{cases}B_3(x,\alpha),&x\geq\xi,\\
B_3(\xi,2x+\alpha-2\xi),&x<\xi,\ 2x+\alpha\geq2\xi,\\
0,&x<\xi,\ 2x+\alpha<2\xi.\end{cases}\label{4.2b}
\end{align}
\end{subequations}
Hence, to solve \eqref{4.1} for any fixed $\xi$, we can limit ourselves
to solving it for $x\geq\xi$, after replacing $B_{1,\xi}$ and $B_{3,\xi}$
with $B_1$ and $B_3$ and $q_\xi$ with $q$.

\subsection{Computation of $\hat{a}_\xi$ and $\hat{b}_\xi$}

As we are considering truncated potentials, we do not need to use
the formulas \eqref{3.14}. Indeed, more simply, we can compute them
by using \eqref{3.17}.
Note that, as a consequence of the decay of $B_1$ and $B_3$, which is
of the same order as for $q$, the above formulas could also be used
numerically
in the case of a non-truncated potential. Let us emphasize that, in
such a way, we are able to compute $\hat{a}_\xi$ and $\hat{b}_\xi$
without resorting to the Zakharov-Shabat system.

\subsection{Computation of $\hat{R}_\xi(\alpha)$}

The computation of $\hat{R}_\xi$ can be carried out by means of three
FFT's. In fact we need two FFT's to compute $a_\xi(\lambda)$ and
$b_\xi(\lambda)$ from $\hat{a}_\xi(\alpha)$ and $\hat{b}_\xi(\alpha)$. We need a
third FFT to compute $\hat{R}_\xi(\alpha)$ from
$R_\xi(\lambda)=-b_\xi(\lambda)/a_\xi(\lambda)$.

\subsection{Computation of the Initial Marchenko kernel}

To compute the Marchenko kernel we still have to evaluate the parameters
of a possible bound state. For this purpose, we resort to an artifice:
we introduce an additional truncated potential $q_\zeta$, with $\zeta>\xi$
so large as not to admit a bound state, which is true if $\|q_\zeta\|_1<\pi/2$.
As a consequence $\Omega_\zeta(\alpha)=\hat{R}_\zeta(\alpha)$, while
$\Omega_\xi(\alpha)=\hat{R}_\xi(\alpha)+\Gamma_\xi e^{-\gamma_\xi\alpha}$, where
$\Gamma_\xi=0$ if and only if there is no bound state.

We now note that for $\alpha\geq2\zeta$ we have
$B_1(\alpha,\,\cdot\,)=B_{1,\xi}(\alpha,\,\cdot\,)
=B_{1,\zeta}(\alpha,\,\cdot\,)$ and
$B_2(\alpha,\,\cdot\,)=B_{2,\xi}(\alpha,\,\cdot\,)
=B_{2,\zeta}(\alpha,\,\cdot\,)$, which implies that both
$\Omega_\xi(\alpha)$ and $\Omega_\zeta(\alpha)$ solve the same
integral system,
so that $\Gamma_\xi
e^{-\gamma_\xi\alpha}=\hat{R}_\zeta(\alpha)-\hat{R}_\xi(\alpha)$.
As a result we can compute numerically
$\hat{R}_\zeta(\alpha)-\hat{R}_\xi(\alpha)$ for $\alpha\geq2\zeta$ and then
apply a least squares technique to estimate the parameters
$\Gamma_\xi$ and $\gamma_\xi$. Hereafter, for the sake of simplicity,
we shall omit the subscript $\xi$.

\subsection{Time evolution of the Marchenko kernel}

Let us assume we have one bound state and the order of the
corresponding zero $\gamma$ of $a(\lambda)$ is 1. In this case
$\Omega(\alpha,t)=\hat{R}(\alpha,t)+\Gamma e^{-\gamma\alpha} e^{-4\mathbf{i}\gamma^2 t}$
where $\hat{R}(\alpha,t)$ is the solution of the initial value
differential problem
\begin{equation*}
\begin{cases}
\frac{\partial\hat{R}}{\partial t}
+4\mathbf{i}\frac{\partial^2\hat{R}}{\partial\alpha^2}=0
, & \alpha\in\mathbb{R},\ t>0 \\[2mm]
\hat{R}(\alpha,0) = \hat{R}(\alpha)
\end{cases}
\end{equation*}

For an effective approximation of $\hat{R}(\alpha,t)$ we distinguish
two situations:
\begin{itemize}
  \item we have to compute $\hat{R}(\alpha,t)$ at a fixed $t=t_0$,
    with $t_0$ not too small ($t_0\geq1$, say);
  \item we have to compute $\hat{R}(\alpha,t)$ for $0<t<t_0$.
\end{itemize}
In the former case we compute $\hat{R}(\alpha,t_0)$ by using \eqref{3.19},
i.e., we compute a discrete convolution. In the latter case the
best strategy consists of discretizing the differential equation with
respect to $\alpha$ and then solving the resulting equation by using
the Discrete Fourier Transform (DFT).

The first step generates the biinfinite system
\begin{equation}\label{4.3}
  \frac{\partial \hat{R}_j}{\partial t}
  =
  -\frac{4\mathbf{i}}{h^2} \big(\hat{R}_{j+1}-2\hat{R}_j+\hat{R}_{j-1}\big),
  \quad
  \hat{R}_j=\hat{R}_j(jh,t),
  \quad
  h>0,\ j\in\mathbb{Z}.
\end{equation}
Using the DFT; i.e., the transform
$$
  \tilde{S}(z) = \mathcal{F}\{S_j\} = \sum_{j=-\infty}^\infty S_j z^j,
  \quad z=e^{\mathbf{i}\theta},
$$
we convert the difference equation \eqref{4.3} into the ordinary
differential equation
$$
  \frac{\partial}{\partial t}\mathcal{F}\{\hat{R}(z,t)\}
  =\frac{8\mathbf{i}}{h^2} (1-\cos(\theta))\mathcal{F}\{\hat{R}(z,t)\}
$$
whose solution is
$$
 \mathcal{F}\{\hat{R}(z,t)\}
= e^{-\frac{4\mathbf{i} t}{h^2}(z+z^{-1}-2)}\mathcal{F}\{\hat{R}(z,0)\},
$$
where $\mathcal{F}\{\hat{R}(z,0)\}$ is known. The computation of
$\mathcal{F}\{\hat{R}(z,t)\}/\mathcal{F}\{\hat{R}(z,0)\}$ can then be done
resorting to a Bessel function and the inverse DFT gives us
$\{\hat{R}(jh,t)\}$.

\subsection{Evolution in time of the initial potential}

As explained in Section 2 the evolution in time of the initial
potential, as the one of the initial energy, is immediately given
(formulas \eqref{2.7} e \eqref{2.8}) by the solution of the Marchenko
system \eqref{2.4}. For fixed values of $t$ and $x$, we discretize
the Marchenko systems at the equispaced points $\{\alpha_i=ih\}$ by
using the composite Simpson rule at the knots $\{\beta_j=jh\}$. In
this way we obtain the linear system
\begin{equation*}
  \left(\begin{array}{cc}
    I  & -H^*D\\
    HD & I
  \end{array}\right)
  \left(\begin{array}{cc}
    \mathbf{b}_1\\
    \mathbf{b}_2
  \end{array}\right)
  =
  \left(\begin{array}{cc}
    \mathbf{0}\\
    -\hat{\mathbf{R}}
  \end{array}\right)
\end{equation*}
where $D=\frac{h}3\mathop{\rm diag}(1,4,2\dots,2,4,1)$, $\mathbf{b}_s=B_s(ih,0,t)$
for $s=1,2$, $H$ is the Hankel matrix defined by $H_{ij}=\Omega(2x+\alpha_i+\beta_j)$ and
$\hat{\mathbf{R}}=\hat{R}(2ih,t)$. This system can be written as
follows
\begin{equation*}
  \begin{cases}
    (D+DHDH^*D)\mathbf{b}_2=-D\mathbf{R},
    \\
    \mathbf{b}_1 = H^* D \mathbf{b}_2,
  \end{cases}
\end{equation*}
where in the first equation the matrix is symmetric positive definite
and the solution is smooth so that it can be easily obtained by the
Conjugate Gradient (CG) method.

As remarked before, the potential $q(x_i,t)$ is given by
$2B_2(x_i,0,t)$ and its energy $\int_{x_i}^{+\infty}|q(y,t)|^2dy$
by $-2B_1(x_i,0,t)$.

\section{Numerical results}\label{sec:5}

To assess the effectiveness of our method, we solved the
NLS equation assuming as the initial potential the one defined in
\eqref{3.15},
where $q(x,0)=-8e^{2x}/(1+4e^{4x})$ is a single lobe potential
\cite{KS03}. More precisely, it is the one soliton potential
considered in Appendix \ref{sec:B}, where $a=p=c=1$. As $\|q\|_1=\pi$, the
truncated potential $q_\xi$ leads to at most one bound state. Since
$\int_{x_0}^{\infty}|q(x,0)|\,dx=\frac\pi2$ (see Appendix \ref{sec:A}
for further details), there is exactly one bound state if and only if
$\xi\leq x_0=-\ln(\sqrt{2})$. The possibility of solving the NLS equation,
by using this initial potential, justifies the choice.  The
analytical solution of the corresponding NLS equation, as well as the
analytical computation of the initial Marchenko kernel $\Omega_\xi$
associated to the above family of truncated potentials $q_\xi$, are
given in Appendix \ref{sec:B}.

In our numerical experiments we considered the two truncated
potentials: $q_{-1}(x,0)$ and $q_1(x,0)$, which lead to one and zero bound
states, respectively. As input data, we used the sampling of the
initial potential in a grid of $2^{10}+1$ equispaced points of the
interval $[-7.5,9]$, roughly speaking is the support of $q$.
We need this choice to make the steps described in Section \ref{sec:4}.

Taking $q_{-1}(x,0)$, the algorithm recognizes automatically the
existence of one bound state and approximates, at a satisfactory
precision, the parameters $\Gamma$ and $\gamma$ identifying the
bound state.  In our experiments we computed the potential in the
rectangle $\mathcal{D}=\{\,(x,t),\text{ s.t. } x\in[-1.90,1.65] \text{
and } t\in[1,1.6]\}$, where the solution is not negligible (note that
in the final step of the algorithm we can choose the subdomain where we
want to compute the potential). The infinite norm of the pointwise
discretization error; i.e., of the difference between the left and
the right hand side of the NLS equation in the grid points of the
domain considered is $\simeq10^{-3}$ if we take $q_{-1}$ as the initial
potential and $\simeq10^{-6}$ in the $q_1$ case (this is mainly due to
the fact that $x_0$ is much further from 1 than from $-1$).

In Figures \ref{fig:-1E} and \ref{fig:1E}, respectively, we depicted
the real and imaginary parts, as well as the modulus, of the
discretization error in $\mathcal{D}$ pertaining to the potentials
$q_{-1}$ and $q_1$, respectively.

\begin{figure}[!ht]
  \centering
  \subfigure[$\Re(\mathbf{i} q_t-q_{xx}-2q|q|^2)$]
  {\includegraphics[width=3.9cm]{fig1a}}
  \quad
  \subfigure[$\Im(\mathbf{i} q_t-q_{xx}-2q|q|^2)$]
  {\includegraphics[width=3.9cm]{fig1b}}
  \quad
  \subfigure[$\big|\mathbf{i} q_t-q_{xx}-2q|q|^2\big|$]
  {\includegraphics[width=3.9cm]{fig1c}}
  \caption{Discretization error determined by $q_{-1}$ (in $\mathcal{D}$)}
  \label{fig:-1E}
\end{figure}
%
\begin{figure}[!ht]
  \centering
  \subfigure[$\Re(\mathbf{i} q_t-q_{xx}-2q|q|^2)$]
  {\includegraphics[width=3.9cm]{fig2a}}
  \quad
  \subfigure[$\Im(\mathbf{i} q_t-q_{xx}-2q|q|^2)$]
  {\includegraphics[width=3.9cm]{fig2b}}
  \quad
  \subfigure[$\big|\mathbf{i} q_t-q_{xx}-2q|q|^2\big|$]
  {\includegraphics[width=3.9cm]{fig2c}}
  \caption{Discretization error determined by $q_{1}$ (in
  $\mathcal{D}$)}
  \label{fig:1E}
\end{figure}

The reflection coefficient $\hat{R}_\xi$ is real (see Appendix
\ref{sec:B}). Our numerical computations give, for both $\xi=1$ and
$\xi=-1$, a componentwise imaginary part which is $\simeq10^{-17}$.
Hence, in Figures \ref{fig:-1O} and \ref{fig:1O} we give the real
part of $\hat{R}_\xi$ and the discretization error in the interval
$[-15,18]$.

\begin{figure}[!ht]
  \centering
  \subfigure[$\Re(\hat{R}_{-1})$]
  {\includegraphics[width=5.9cm]{fig3a}}
  \quad
  \subfigure[$\big|\hat{R}_{-1}-\hat{R}_{-1}^{true}\big|$]
  {\includegraphics[width=5.9cm]{fig3b}}
  \caption{$\Re\big(\hat{R}_{-1}(x,0)\big)$ and the associated pointwise error}
  \label{fig:-1O}
\end{figure}

\begin{figure}[!ht]
  \centering
  \subfigure[$\Re(\hat{R}_{1})$]
  {\includegraphics[width=5.9cm]{fig4a}}
  \quad
  \subfigure[$\big|\hat{R}_{1}-\hat{R}_{1}^{true}\big|$]
  {\includegraphics[width=5.9cm]{fig4b}}
  \caption{$\Re\big(\hat{R}_{1}(x,0)\big)$ and the associated pointwise error}
  \label{fig:1O}
\end{figure}

Finally in Figures \ref{fig:-1T} and \ref{fig:1T} we visualize the
time evolution of $q_{-1}$ and $q_1$, respectively (in $\mathcal{D}$).

\begin{figure}[!ht]
  \centering
  \subfigure[$\Re(q_{-1})$]
   {\includegraphics[width=5.9cm]{fig5a}}
  \quad
  \subfigure[$\Im(q_{-1})$]
   {\includegraphics[width=5.9cm]{fig5b}}
  \caption{Real and imaginary parts of $q_{-1}(x,t)$}
  \label{fig:-1T}
\end{figure}

\begin{figure}[!ht]
  \centering
  \subfigure[$\Re(q_1)$]
  {\includegraphics[width=5.9cm]{fig6a}}
  \quad
  \subfigure[$\Im(q_1)$]
  {\includegraphics[width=5.9cm]{fig6b}}
  \caption{Real and imaginary parts of $q_{1}(x,t)$}
  \label{fig:1T}
\end{figure}

\section{Conclusions}

In this paper we give a new method to solve the NLS equation \eqref{1.1}.
Though it follows the path of the IST, it works without solving the
Za\-kha\-rov-Sha\-bat system. From the numerical point of view it
uses the structure of the kernels in the Marchenko integral equations
optimally. The effectiveness of our method has been ascertained under
the hypothesis that the initial potential has at most one bound
state. It is an open problem to develop a method to treat the
general case of $n$ bound states. A second open problem, of great relevance
to telecommunications, concerns the extension of our method to the
solution of NLS equation with a damping factor in the energy term.

\section{Number of Bound States}\label{sec:A}

In this appendix we summarize the results in \cite{KS03} concerning the number
of bound states:
\begin{itemize}
\item[(1)] There are no bound states if $\|q\|_1\le(\pi/2)$.

\item[(2)] There are no bound states if $q$ is real and
$$
\|q\|_1\le\max\left(\|\max(q,0)\|_1,\|\max(-q,0)\|_1\right)
\le\frac{\pi}{2}.
$$

\item[(3)] If $q$ is a real potential of constant sign such that
\begin{equation}\label{A.1}
(2N-1)\frac{\pi}{2}<\|q\|_1\le(2N+1)\frac{\pi}{2},
\end{equation}
then there are at least $N$ positive imaginary bound state spectral parameters
and no bound state spectral parameter in $\mathbb{C}^+\setminus(\mathbf{i}\mathbb{R}^+)$ of maximal
imaginary part among the bound state spectral parameters in $\mathbb{C}^+$.

\item[(4)] If $q$ is a single lobe potential [i.e., if $q$ is real and has
constant sign, and $|q|$ is nonincreasing on $\mathbb{R}^+$ and nondecreasing on
$\mathbb{R}^-$], then the number of positive imaginary bound state spectral parameters
is exactly $N$ and there are no bound state spectral parameters off the
imaginary axis.

\end{itemize}
Thus if both $\int_x^\infty dy\,\max(q(y,0),0)<(\pi/2)$ and
$\int_x^\infty dy\,\max(-q(y,0),0)<(\pi/2)$, then $p(x)=0$. Further, if $q$ is
a single lobe potential satisfying
$$
(2N-1)\frac{\pi}{2}<\int_x^\infty dy\,|q(y,0)|\le(2N+1)\frac{\pi}{2}
$$
for some positive integer $N$, then $p(x)=N$.

\section{Truncated One Soliton Potentials}\label{sec:B}

In this appendix we give the Jost solutions and scattering coefficients of the
truncated one-soliton potential.

Consider the initial potential
$$
q(x,0)=\frac{-2ce^{-2ax}}{1+\frac{|c|^2}{4p^2}e^{-4px}},\quad x\in\mathbb{R},
$$
where $a=p+\mathbf{i} q$ with $p>0$ and $0\neq c\in\mathbb{C}$, which stands
for one soliton potential. The application of direct and inverse
scattering transform leads to the following data:
\begin{gather*}
\|q\|_1=\pi, \quad\Omega(\alpha)=ce^{-a\alpha},\\
B_1(x,\alpha)=\frac{-\frac{|c|^2}{2p}e^{-4px}e^{-\overline{a}\alpha}}
{1+(|c|^2/4p^2)e^{-4px}},\quad
B_2(x,\alpha)=\frac{-ce^{-2ax}e^{-a\alpha}}{1+(|c|^2/4p^2)e^{-4px}},\\
\Psi_1(\lambda,x)=e^{\mathbf{i}\lambda x}\frac{\lambda-\mathbf{i}\gamma(x)}{\lambda+\mathbf{i}\,\overline{a}},\quad
\Psi_2(\lambda,x)=\frac{\mathbf{i} ce^{-2ax}e^{-\mathbf{i}\lambda x}}
{(\lambda-ia)[1+(|c|^2/4p^2)e^{-4px}]},
\end{gather*}
where
$$
\gamma(x)=p\frac{-1+\frac{|c|^2}{4p^2}e^{-4px}}{1+(|c|^2/4p^2)e^{-4px}}+\mathbf{i} q.
$$
Moreover, as it is easy to check, the NLS solution is given by
$$
q(x,t)=\frac{-2c\,e^{-2p(x-4qt)}e^{-2\mathbf{i} q(x-2qt,0)}e^{-4\mathbf{i} tp^2}}
{1+(|c|^2/4p^2)e^{-4p(x-4qt)}}.
$$

Now use subscripts $\xi$ to denote quantities pertaining to the truncated
potential $q_\xi$. Then the functions which characterize the
reflection coefficient $R(\lambda)$ are:
\[
a_\xi(\lambda)=e^{-\mathbf{i}\lambda\xi}\Psi_1(\lambda,\xi)
=\frac{\lambda-\mathbf{i}\gamma(\xi)}{\lambda+\mathbf{i}\overline{a}}=1-\int_0^\infty d\alpha\,
[\gamma(\xi)+\overline{a}]e^{\mathbf{i}\lambda\alpha}e^{-\overline{a}\alpha},
\]
and
\begin{align*}
b_\xi(\lambda)&=e^{-\mathbf{i}\lambda\xi}\Psi_2(\lambda,\xi)\nonumber\\
&=\frac{\mathbf{i} ce^{-2a\xi}e^{-2\mathbf{i}\lambda\xi}}
{(\lambda-\mathbf{i} a)[1+(|c|^2/4p^2)e^{-4p\xi}]}\\
&=\frac{-c}{1+(|c|^2/4p^2)
e^{-4p\xi}}\int_{2\xi}^\infty d\alpha\,e^{-\mathbf{i}\lambda\alpha}e^{-a\alpha}.
\end{align*}
Note that $\mathbf{i}\gamma(\xi)\in\mathbb{C}^+$ if and only if
$\xi<x_0:=[\ln(|c|/2p)]/2p$, which implies that
we have a bound state if and only if $\xi<x_0$. Therefore,
$$
R_\xi(\lambda)
=-\frac{b_\xi(\lambda)}{a_\xi(\lambda)}=\frac{-\mathbf{i} ce^{-2a\xi}e^{-2\mathbf{i}\lambda\xi}}
{1+(|c|^2/4p^2)e^{-4p\xi}}\,\frac{\lambda+\mathbf{i}\overline{a}}{\lambda-\mathbf{i} a}\,\frac{1}
{\lambda-\mathbf{i}\gamma(\xi)}
$$
and
$$
 T_\xi(\lambda)=\frac{\lambda+\mathbf{i}\overline{a}}{\lambda-\mathbf{i}\gamma(\xi)}.$$
Let us now apply the inversion formula
$$
\hat{R}_\xi(\alpha)=\frac{ce^{-2a\xi}}{1+(|c|^2/4p^2)e^{-4p\xi}}\,
\frac{1}{2\pi\mathbf{i}}\int_{-\infty}^\infty d\lambda\,e^{\mathbf{i}\lambda(\alpha-2\xi)}
\frac{\lambda+\mathbf{i}\overline{a}}{\lambda-\mathbf{i} a}\,\frac{1}{\lambda-\mathbf{i}\gamma(\xi)}.
$$
 For $\xi>x_0$ we get
$$
\Omega_\xi(\alpha)=\hat{R}_\xi(\alpha)=\begin{cases}\frac{2pc}{[a-\gamma(\xi)]
[1+(|c|^2/4p^2)e^{-4p\xi}]}e^{-a\alpha}=ce^{-a\alpha},&\alpha>2\xi,\\[3pt]
\frac{-ce^{2\xi[\gamma(\xi)-a]}}{1+(|c|^2/4p^2)e^{-4p\xi}}\,
\frac{\gamma(\xi)+\overline{a}}{\gamma(\xi)-a}\,e^{-\gamma(\xi)\alpha},&\alpha<2\xi.
\end{cases}
$$
On the other hand, for $\xi<x_0$ we have $N=1$, $n_1=1$, $\lambda_1=\mathbf{i}\gamma(\xi)$,
and
$$
\Gamma(\xi)=-\frac{c}{1+(|c|^2/4p^2)e^{-4p\xi}}\,
\frac{\gamma(\xi)+\overline{a}}{\gamma(\xi)-a}e^{2\xi[\gamma(\xi)-a]},
$$
whence
$$
\hat{R}_\xi(\alpha)=\begin{cases}ce^{-a\alpha}
-\Gamma(\xi)e^{-\gamma(\xi)\alpha},&\alpha>2\xi,\\ 0,&\alpha<2\xi.\end{cases}
$$
As a result, for $\xi<x_0$ we have
$$
\Omega_\xi(\alpha)=\hat{R}_\xi(\alpha)+\Gamma_{10}e^{-\gamma(\xi)\alpha}=\begin{cases}
ce^{-a\alpha},&\alpha>2\xi,\\ \Gamma_{10}e^{-\gamma(\xi)\alpha},&\alpha<2\xi.\end{cases}
$$
Thus in either case ($\xi<x_0$ and $\xi>x_0$) we have found the same
expression for $\Omega_\xi(\alpha)$, both resorting to the Zakharov-Shabat
system and ignoring it.

Let us compute now
\begin{align*}
\int_\xi^\infty dx\,|q(x,0)|&=2|c|\int_\xi^\infty dx\,\frac{e^{-2px}}
{1+(|c|^2/4p^2)e^{-4px}}\\ &\hskip-25pt\stackrel{z=(|c|/2p)e^{-2px}}{=}
2\int_0^{(|c|/2p)e^{-2p\xi}}\frac{dz}{1+z^2}
=2\arctan\Big(\frac{|c|}{2p}e^{-2p\xi}\Big).
\end{align*}
Thus $\int_\xi^\infty dx\,|q(x,0)|=(\pi/2)$ if and only if $\xi=x_0$
which implies that $q_\xi(x)$ has one bound state if $\xi\leq x_0$ and
no bound state if $\xi>x_0$. Also,
$$
\int_{-\infty}^\infty dx\,|q(x,0)|^2=\big[\frac{4p}{1+(|c|^2/4p^2)
e^{-4px}}\big]_{x=-\infty}^\infty=4p.
$$

\begin{thebibliography}{WW}

   \bibitem{AC}
     M. J. Ablowitz and P. A. Clarkson,
     \emph{Solitons, Nonlinear Evolution Equations and Inverse Scattering},
     Cambridge Univ. Press, Cambridge, 1991.

   \bibitem{AKNS}
     M. J. Ablowitz, D. J. Kaup, A. C. Newell, and H. Segur,
     \emph{The inverse scattering transform-Fourier analysis for nonlinear
     problems},
     Stud. Appl. Math. \textbf{53}, 249--315 (1974).

   \bibitem{AL76}
     M. J. Ablowitz and J. F. Ladik,
     \emph{A nonlinear difference scheme and inverse scattering},
     Studies in Appl. Math. \textbf{55}, 213--229 (1976).

   \bibitem{AL77}
     M. J. Ablowitz and J. F. Ladik,
     \emph{On the solution of a class of nonlinear partial difference equations},
     Studies in Appl. Math. \textbf{57}, 1--12 (1976/77).

   \bibitem{APT}
     M. J. Ablowitz, B. Prinari, and A. D. Trubatch,
     \emph{Discrete and Continuous Nonlinear Schr\"o\-din\-ger Systems},
     Cambridge Univ. Press, Cambridge, 2004.

   \bibitem{AS81}
     M. J. Ablowitz and H. Segur,
     \emph{Solitons and the Inverse Scattering Transform},
     SIAM, Philadelphia, 1981.

   \bibitem{AS64}
     M. Abramowitz and I. A. Stegun,
     \emph{Handbook of Mathematical Physics},
     Dover Publ., New York, 1964.

   \bibitem{Ag01} G. P. Agrawal,
     \emph{Nonlinear Fiber Optics},
     3rd ed., Academic Press, New York, 2001.

   \bibitem{ADV07}
     T. Aktosun, F. Demontis, and C. van der Mee,
     \emph{Exact solutions to the focusing nonlinear Schr\"o\-din\-ger equation},
     Inverse Problems \textbf{23}, 2171--2195 (2007).

   \bibitem{AKV00}
     T. Aktosun, M. Klaus, and C. van der Mee,
     \emph{Direct and inverse scattering for selfadjoint Hamiltonian system
     on the line},
     Integral Equations and Operator Theory \textbf{38}, 129--171 (2000).

   \bibitem{B}
     Th. Busse,
     \emph{Generalized Inverse Scattering Transform for the Nonlinear
     Schr\"o\-din\-ger Equation},
     Ph.D. Thesis, University of Texas at Arlington, 2008.

   \bibitem{C03}
     Th. Cazenave,
     \emph{Semilinear Schr\"o\-din\-ger Equations},
     Courant LNM \textbf{10}, Amer. Math. Soc., Providence, RI, 2003.

   \bibitem{C73}
     J. B. Conway,
     \emph{Functions of One Complex Variable},
     Graduate Texts in Mathematics \textbf{11}, Springer, Berlin, 1973.

   \bibitem{D07}
     F. Demontis, \emph{Direct and Inverse Scattering of the Matrix
     Zakharov-Shabat System},
     Ph.D. Thesis, University of Cagliari, Italy, 2007.

   \bibitem{DEGM}
     R. K. Dodd, J. C. Eilbeck, J. D. Gibbon, and H. C. Morris,
     \emph{Solitons and Nonlinear Wave Equations},
     Academic Press, London, 1982.

   \bibitem{EH}
     W. Eckhaus and A. van Harten,
     \emph{The Inverse Scattering Transformation and the Theory of
     Solitons},
     North-Holland, Amsterdam, 1981.

   \bibitem{F96}
     B. Fornberg,
     \emph{A Practical Guide to Pseudospectral Methods},
     Cambridge University Press, Cambridge, 1996.

   \bibitem{FW73}
     B. Fornberg and G. B. Whitham,
     \emph{A Numerical and theoretical study of certain nonlinear wave
     phenomena},
     Philos. Trans. Roy. Soc. London, Series A, \textbf{289}(1361), 373--404
     (1973).

   \bibitem{GGKM67}
     C. S. Gardner, J. M. Greene, M. D. Kruskal, and R. M. Miura,
     \emph{Method for solving the Kor\-te\-weg-de Vries equation},
     Phys. Rev. Lett. \textbf{19}, 1095--1097 (1967).

   \bibitem{GGKM74}
     C. S. Gardner, J.M. Greene, M. D. Kruskal, and R. M. Miura,
     \emph{The Kor\-te\-weg-de Vries equation and
     generalizations. VI. Methods for exact solution},
     Commun. Pure Appl. Math. \textbf{27}, 97--133 (1974).

   \bibitem{HaT73}
     R. H. Hardin and F. D. Tappert,
     \emph{Applications of the split-step Fourier method to the
     numerical solution of non-linear and variable-coefficient wave
     equations},
     SIAM-SIGNUM Fall Meeting, Austin, October, 1972, SIAM
     Rev. Chronicle \textbf{15}, 423 (1973).

   \bibitem{HM02}
     A. Hasegawa and M. Matsumoto,
     \emph{Optical Solitons in Fibers},
     3rd ed., Springer, Berlin, 2002.

   \bibitem{HT73}
     A. Hasegawa and F. Tappert,
     \emph{Transmission of stationary nonlinear optical pulses in
     dispersive dielectric fibers. I. Anomalous dispersion},
     Appl. Phys. Lett. \textbf{23}, 142--144 (1973).

   \bibitem{KS03}
     M. Klaus and J. K. Shaw,
     \emph{On the eigenvalues of Zakharov-Shabat systems},
     SIAM J. Math. Anal. \textbf{34}(4), 759--773 (2003).

   \bibitem{LYRF}
     B. M. Lake, H.C. Yuen, H. Rungaldier, and W. E. Ferguson,
     \emph{Nonlinear deep-water waves; theory and experiment. Part} 2.
     \emph{Evolution of a continuous wave train},
     J. Fluid Mech. \textbf{83}, 49--74 (1977).

   \bibitem{M86}
     V. A. Marchenko,
     \emph{Sturm-Liouville Operators and Applications},
     Birk\-h\"au\-ser OT\textbf{22}, Basel and Boston, 1986; also:
     Naukova Dumka, Kiev, 1977 [Russian].

   \bibitem{NMPZ}
     S. Novikov, S. V. Manakov, L. P. Pitaevskii, and V. E. Zakharov,
     \emph{Theory of Solitons: The Inverse Scattering Method},
     Consultants Bureau, New York, 1984.

   \bibitem{R97}
     M. P. Robinson,
     \emph{The solution of nonlinear Schr\"o\-din\-ger equations using
     orthogonal spline collocation},
     Comput. Math. Appl. \textbf{33}, 39--57 (1997);
     also: Comput. Math. Appl. \textbf{35}, 151 (1998).

   \bibitem{RF94}
     M. P. Robinson and G. Fairweather,
     \emph{Orthogonal spline collocation methods for
     Schr\"o\-din\-ger-type equations in one space
     variable},
     Numer. Math. \textbf{68}, 355--376 (1994).

   \bibitem{RFH93}
     M. P. Robinson, G. Fairweather, and B. M. Herbst,
     \emph{On the numerical solution of the cubic Schr\"o\-din\-ger
     equation in one space variable},
     J. Comput. Phys. \textbf{104}, 277--284 (1993).

   \bibitem{SCM}
     A. C. Scott, F. Y. F. Chu, and D. W. McLaughlin,
     \emph{The soliton: A new concept in applied science},
     Proc. IEEE \textbf{61}, 1443--1483 (1973).

   \bibitem{Sh04}
     J. K. Shaw,
     \emph{Mathematical Principles of Optical Fiber Communications},
     CBMS-NSF Regional Conference Series \textbf{76}, SIAM, Philadelphia, 2004.

   \bibitem{TA84a}
     Th. R. Taha and M. J. Ablowitz,
     \emph{Analytical and numerical aspects of certain nonlinear
     evolution equations}. I.~\emph{Analytical},
     J. Comput. Phys. \textbf{55}, 192--202 (1984).

   \bibitem{TA84b}
     Th. R. Taha and M. J. Ablowitz,
     \emph{Analytical and numerical aspects of certain nonlinear
     evolution equations}. II.~\emph{Numerical, nonlinear
     Schr\"o\-din\-ger equation},
     J. Comput. Phys. \textbf{55}, 203--230 (1984).

   \bibitem{V04}
     C. van der Mee \emph{Direct and inverse scattering for
     skewselfadjoint Hamiltonian systems}.
     In: J.A. Ball, J.W. Helton, M. Klaus, and L. Rodman (eds.),
     \emph{Current Trends in Operator Theory and its Applications},
     Birk\-h\"au\-ser OT \textbf{149}, Basel and Boston, 2004, pp. 407--439.

   \bibitem{VAP05}
     J. Villarroel, M. J. Ablowitz, and B. Prinari,
     \emph{Solvability of the direct and inverse problems for the
     nonlinear Schr\"o\-din\-ger equation},
     Acta App. Math. \textbf{87}, 245--280 (2005).

   \bibitem{WH86}
     J. A. C. Weideman and B. M. Herbst,
     \emph{Split-step methods for the solution of the nonlinear
     Schr\"o\-din\-ger equation},
     SIAM J. Numer. Anal. \textbf{23}, 485--507 (1986).

   \bibitem{W74}
     G. B. Whitham,
     \emph{Linear and Nonlinear Waves}, John Wiley, New York, 1974.

   \bibitem{XT03}
     X. Xu and T. Taha,
     \emph{Parallel split-step Fourier methods for nonlinear
     Schr\"o\-dinger-type equations},
     J. Math. Model. Algorithms 2, (2003),  \textbf{3}, 185--201.

   \bibitem{YL75}
     H. C. Yuen and B. M. Lake,
     \emph{Nonlinear deepwater waves}: \emph{Theory and experiment},
     Phys. Fluids \textbf{18}, 756--760 (1975).

   \bibitem{Z68}
     V. E. Zakharov,
     \emph{Stability of periodic waves of finite amplitude on the
     surface of a deep fluid},
     J. Appl. Mech. Tech. Phys. \textbf{4}, 190--194 (1968).

   \bibitem{ZS}
     V. E. Zakharov and A. B. Shabat,
     \emph{Exact theory of two-dimensional self-focusing and one
     dimensional self-modulation of waves in nonlinear media},
     Sov. Phys. JETP \textbf{34}, 62--69 (1972).

\end{thebibliography}

\end{document}