An Approach to Solving Direct and Inverse Scattering Problems for Non-Selfadjoint Schrödinger Operators on a Half-Line

In this paper, an approach to solving direct and inverse scattering problems on the half-line for a one-dimensional Schrödinger equation with a complex-valued potential that is exponentially decreasing at infinity is developed. It is based on a power series representation of the Jost solution in a unit disk of a complex variable related to the spectral parameter by a Möbius transformation. This representation leads to an efficient method of solving the corresponding direct scattering problem for a given potential, while the solution to the inverse problem is reduced to the computation of the first coefficient of the power series from a system of linear algebraic equations. The approach to solving these direct and inverse scattering problems is illustrated by several explicit examples and numerical testing.

Studying a Zakharov-Shabat system, even with a real-valued potential, naturally leads to a couple of equations of the form (1) with complex-valued potentials; see [9]. Indeed, consider the Zakharov-Shabat system where ρ is a complex spectral parameter and u(x) is a real-valued potential.
This solution is called the Jost solution of (1). It admits the Levin integral representation [12] (see also [18,23,24]) where for every fixed x, the kernel A(x, t) belongs to L 2 (x, ∞). In [25] (see also [26]) a Fourier-Laguerre series representation for A(x, t) was proposed in the form A(x, t) = ∞ ∑ n=0 a n (x)L n (t − x)e where L n (τ) stands for the Laguerre polynomial of order n. A recurrent integration procedure was developed in [27] to calculate the coefficients a n (x). The substitution of (7) into (6) was found to lead to a series representation for the Jost solution [25,26] e(ρ, x) = e iρx 1 + (z + 1) ∞ ∑ n=0 (−1) n z n a n (x) , x ≥ 0, ρ ∈ C + where In the present work, we consider the direct and inverse scattering problems for (1) subject to the homogeneous Dirichlet condition however, the approach developed here is also applicable in the case of other boundary conditions, such as y (0) − hy(0) = 0 with h ∈ C.
The problem (1) and (10) under Condition (2) possesses a continuous spectrum coinciding with the positive semi-axis λ > 0, and may have a point spectrum that coincides with the squares of the non-real roots of the Jost function e(ρ) := e(ρ, 0), if such roots exist. Let us denote them as ρ 1 , . . . , ρ α . Their multiplicity may be greater than one. In this case, instead of norming constants associated to the eigenvalues, the corresponding normalization polynomials X k (x) naturally arise (see Section 3.3 below).
As a component of the scattering data for (1), the scattering function is considered in the strip |Im(ρ)| < ε 0 where ε 0 is sufficiently small (see Section 3.2 below). The direct scattering problem for (1) and (10) consists of obtaining the set of the scattering data {ρ k , m k , X k (x)} α k=1 , s(ρ) .
The overall approach developed in the present work to solve this problem is based on the representation (8). Indeed, the calculation of {ρ k } α k=1 is easily realizable with the aid of the argument principle theorem applied to find zeros of (8) in the unit disc. To the best of our knowledge, there has been no practical way of calculating the normalization polynomials. We propose a simple procedure for computing their coefficients by solving a finite system of linear algebraic equations. For this, an auxiliary result for the derivatives ∂ m ∂z m e(ρ(z), x) is obtained.
The calculation of the scattering function s(p) requires an analytic extension of the Jost function e(ρ) obtained from (8), onto the strip −ε 0 < Im(ρ) < 0. We explore different possibilities for such an extension, including the Padé approximants (see [28,29]) and the power series analytic continuation [30] (p. 150), [31]. This results in an efficient numerical method for solving the direct scattering problem.
The inverse scattering problem consists of recovering the potential q(x) from the set of the scattering data. A general theory of this inverse problem can be found in [12,13,20] (p. 353), [24,[32][33][34][35]. Here, we use the representation (7) for the numerical solution of the problem, thus extending the approach developed in [25,26,[36][37][38] to the non-selfadjoint situation. The inverse Sturm-Liouville problem is reduced to an infinite system of linear algebraic equations. The potential q(x) is recovered from the first component of the solution vector, which coincides with a 0 (x) in (7).
The reduction to the infinite system of linear algebraic equations is based on the substitution of the series representation (7) for the kernel A(x, t) into the Gelfand-Levitan equation (see [39]), where the function f can be computed from the set of scattering data (11): To approximate the complex-valued function a 0 (x), we consider the truncated system of linear algebraic equations, for which the existence, uniqueness and stability of the solution is proved.
Finally, we illustrate the proposed approach by numerical calculations performed in Matlab2021a.
We discuss the details of the numerical implementation of the method: its convergence, stability and accuracy. In a couple of examples, we show the "in-out" performance of the approach, i.e., we solve the direct problem numerically and use the results of our computation as the input data to solve the inverse problem.
The approach based on the representations (7) and (8) leads to efficient numerical methods for solving both direct and inverse scattering problems.
In Section 2, we recall the series representations for the kernel A(x, t) and for the Jost solution, then prove additional results related to these representations. In Section 3, we recall the set of scattering data and put forward an algorithm for solving the direct scattering problem. Additionally, we present analytical examples. In Section 4, the approach for solving the inverse scattering problem is developed. Analytical examples from Section 2 are considered in order to illustrate the approach. In Section 5, we discuss the numerical implementation of the algorithms proposed for solving the direct and inverse scattering problems. Section 6 contains some concluding remarks.

Series Representations for the Transmutation Operator Kernel and Jost Solution
Consider the one-dimensional Schrödinger equation on the half-line (1) where λ = ρ 2 ∈ C is the spectral parameter. The potential q(x) is a complex-valued function satisfying Condition (2) for some ε > 0.

Remark 2. Under Condition
provided the existence of these derivatives; see [21].
Additionally, the kernel A(x, t) has first continuous derivatives that satisfy the inequalities [12] and the equality [12] (p. 328) As was pointed out in [25], since A(x, ·) ∈ L 2 (x, ∞), the function belongs to L 2 ([0, ∞); e −t ) and hence admits the series representation a(x, t) = ∞ ∑ n=0 a n (x)L n (t), (22) where L n (t) stands for the Laguerre polynomial of order n and a n (x) are complex-valued functions such that {a n (x)} ∞ n=0 ∈ l 2 for any x ≥ 0. For all x ≥ 0, the series (22) converges in the norm of L 2 ([0, ∞); e −t ). Thus, and This series representation was obtained in [25] for real-valued q(x). However, (23) remains true in the non-selfadjoint case as well.

Proposition 2.
For any fixed x ≥ 0, the series converges pointwise.
Proof. We use [40] (Theorem 6.5), and thus need to verify that the following assertions are true.

The integrals
exist.
To prove the first assertion, it is enough to consider estimate (15). Indeed, Thus, The second assertion follows from the inclusion A(x, ·) ∈ C 1 (x, ∞). The existence of the first integral in (26) follows from the continuity of a(x, ·). Finally, for the second integral we have and thus, from the proof of the first assertion, we obtain Now, the application of Theorem 6.5 from [40] completes the proof.
The series (8) is convergent in the open unit disk of the complex z-plane, D := {z ∈ C : |z| < 1}, and for every x, the function e(ρ, x)e −iρx belongs to the Hardy space H 2 (D) as a function of z [26]. Proposition 3. Let q(x)(1 + x) ∈ L(0, ∞). Then, the kernel A(x, t) admits the representation (23), where for any x fixed the series converges in the norm of L 2 (x, ∞), and the complex-valued coefficients a n (x) satisfy the system of equations −l[a n ] − a n = −l[a n−1 ] + a n−1 , n = 1, 2, . . . , (29) as well as the inequality Proof. The proof of (28) and (29) from [26] (Theorem 10.1, p. 66) given for the case of a real-valued q remains valid in this more general situation as well.

Remark 3.
Under the assumption that functions a (ν) (x, t) are absolutely continuous with respect to t in [0, ∞) for ν = 0, 1, 2, the convergence of the power series in (8) for z ∈ D can be proved with the aid of a result from [42], which states that and Moreover, .
To ensure Condition (33) for j = 0, notice that from (16) we have For j = 1, Condition (33) holds due to (19). However, the fulfillment of (33) for j = 2 as well as that of (34) requires the additional regularity of q(x), ensuring the possibility of the differentiation of the integral equation for the kernel at least three times [12] (p. 296).

Remark 4.
Denote In [27], the following statements were proved in the case of a real-valued potential.
1. If Im ρ > 0, then These results remain valid in the case of a complex-valued potential. Moreover, under the assumptions of Remark 3, we obtain the inequality .

Remark 5.
The substitution of ρ = i 2 into (8) leads to the equality a 0 (x) = e i 2 , x e x/2 − 1. Moreover, note that we have By ω(ρ, x), we denote the solution of (1), satisfying the initial conditions We also need the solution

Definition 1.
We call the roots of e(ρ) that lie in C + \ {0} the singular numbers of the problem (1) and (10).
If they exist, their number is finite. Let us denote the non-real singular numbers by ρ 1 , . . . , ρ α . The numbers λ k = ρ 2 k constitute the point spectrum of the problem, and the multiplicities of the zeros ρ k (k = 1, . . . , α) are called the multiplicities of the singular numbers and denoted by m k , respectively.
Thus, we are interested in the zeros z k of the Jost function e(ρ) = 1 + (z + 1) ∞ ∑ n=0 (−1) n z n a n (0) (41) to obtain the eigenvalues from λ k = − z k −1 For an estimate of the number of the eigenvalues, we refer to [43].
A function satisfying properties 2-5 is said to be of S-type in the strip |Im ρ| < ε 0 . The following examples illustrate some of the above definitions.
Example 1 ([44,45]). Consider the potential (2). With the aid of Wolfram Mathematica v.12 the Jost solution can be obtained in a closed form, where J ν (z) stands for the Bessel function of the first kind of order ν. Hence, and the eigenvalues are the squares of the values ρ ∈ C + such that From here, we obtain the only singular number The scattering function has the form It is well-defined in the domain and is an S-type function in the strip |Im(ρ)| < 1 2 .
The corresponding Jost solution e 2 (ρ, x) is obtained from the Jost solution of a Zakharov-Shabat system (see [46]) with the potential u(x), Thus, the Jost function is It has one root, ρ * = −1, which corresponds to the spectral singularity λ * = ρ 2 * = 1. The scattering function is given by which is an S-type function in the strip |Im(ρ)| < 1 2 .

Example 3 ([21]
). Consider the potentials of the form satisfying Condition (2) for 0 < ε < 2a. The Jost solution has the form from which the Jost function is obtained The square of this ρ represents the discrete spectrum of the problem. The potential (49) is complex-valued when b is not purely imaginary. The scattering function has the form which is an S-type function in |Im(ρ)| < min{a, Im(−ia tanh(b))} in the case of a complexvalued potential. In the case of a real-valued potential, s(ρ) is an S-type function in |Im(ρ)| < a. To present an explicit example, we fix a = 1 and b = −1 − i in (49). Then, (2) and the Jost solution is Thus, the Jost function has the form and one eigenvalue exists: λ = − tanh 2 (1 + i). The scattering function is an S-type function in the strip |Im(ρ)| < 1.

Normalization Polynomials
The normalization polynomial X k (x) of degree m k − 1, associated with the eigenvalue ρ 2 k (m k is the algebraic multiplicity of ρ k as zero of e(ρ)), defined by the equation [18] i Res(Ω(ρ, x); ρ k ) = e iρ k x X k (x) where Ω(ρ, x) is defined by (40). Using the series representation (23) of the kernel A(x, t), we can obtain a method to compute the coefficients of X k (x).

Remark 6.
Note that the series (8) can be written as in terms of the Jacobi polynomials P (α,β) n (τ).
Let us write Equation (51) in terms of the Jost solution and Jacobi polynomials, as follows.
. . , α be an eigenvalue of problem (1) and (10) and m k be its multiplicity. For the normalization polynomial X k (x), the equality holds Proof. The substitution of (23) into (51) yields Here, we change the order of summation and integration due to Parseval's identity [47] (p. 16) and additionally use the equality Thus, The last integral can be explicitly evaluated [48] (Formula 7.414 (7)) where F(a, b, c; z) stands for the hypergeometric function [49] (p. 56). Thus, we have the equation and due to Remark 6, we obtain (53).
Hereinafter C n k = n k denotes the binomial coefficient.
Proof. We use the identity [50] where F z means the derivative with respect to z, and j, n are integers. Let us prove the lemma by induction. For m = 1, from (52), we have ∂ ∂z (e(ρ, x)) = e iρx x (z + 1) 2 1 + (z + 1) The application of (55) gives Consider Formula (54) as the induction hypothesis for m = k. The idea is to prove the equation and the equality Then, noting that the second terms on the left-hand side of (57) and (58) coincide up to the sign, the desired result is obtained by summing up both equations.
The proof of Equations (57) and (58) is presented in Appendix A, which completes the proof of the Lemma.
As long as there is no possible misunderstanding, we consider a fixed ρ = ρ k with a multiplicity m = m k and the corresponding normalization polynomial X(x) = X k (x). Thus, the index k is omitted along the following two statements.
satisfy the equation Proof. Comparing (53) with (60) we see that, in fact, we need to prove the equality Note that Then, upon comparison of (61) with (62), it can be observed that proving (61) is equivalent to proving the equality for some natural number s ≤ m − 1. Thus, we are going to prove (64). The substitution of the term with the derivative in (63) by Formula (54) for m = r + 1 is enough to obtain (64) as follows This completes the proof of the Lemma.
Equation (60) provides us with a simple method for computing the coefficients b j in (59), and consequently for calculating the normalization polynomials.
to a complex singular number ρ k satisfy the system of linear algebraic equations where A is an m × m k matrix with entries defined by Here, x j ≥ 0 are distinct points, j = 1, . . . , m (m ≥ m k ). B is an m k vector with its entries being the normalization polynomial coefficients B n = b n−1 , n = 1, . . . , m k , and D is an m vector defined by Proof. The proof consists of observing that each row in (65) is just Formula (60) corresponding to a point x j . The number of rows must be at least m k ; otherwise, the system (65) is underdetermined.
Thus, the coefficients of the normalization polynomial are obtained from the system (65).
is called the scattering data set of problem (1) and (10).
Here, ρ k are the non-real singular numbers, m k their multiplicities, X k (x) the corresponding normalization polynomials, and s(ρ) is the scattering function (S-type function in the strip |Im ρ| < ε 0 ).
In order to recall a result on the characterization of the scattering data, we need the following definition [39]. Definition 3. Let s(ρ) be an S-type function in the strip |Im ρ| < ε 0 and let L be a curve lying in the strip and running from −∞ to +∞, such that all roots (poles) of s(ρ) are situated above (below) L. The increment divided by 2π of a continuous branch of Arg s(ρ), when ρ runs along L from −∞ to +∞, is called the index of s(ρ) and denoted by Ind s.
Let us assume that a set J as in Definition 2 is given. A necessary and sufficient condition (obtained in [18]) to ensure that this set represents the scattering data for a problem (1) and (10) with Condition (2) is the following relation In the case when m k = 1, the notion of the Birkhoff solution is useful for computing the corresponding norming constants. (1) (see [24] (p. 113)), i.e., a solution satisfying the asymptotic relation

Remark 7. Let E(ρ, x) denote the Birkhoff solution of Equation
is also a Birkhoff solution of (1) for any constant c ∈ C. Note that for ρ = ρ k (a singular number of the problem), the values of all Birkhoff solutions at the origin coincide. We have E( Note that ρ k is a pole of Ω(ρ, x) in the upper half-plane of the complex variable ρ if and only if it is a root of the Jost function e(ρ) (see (40)). Thus, in case of a simple pole ρ k in Equation (67), the residue can be computed as follows whereė(ρ) := d dρ e(ρ), and the corresponding normalization polynomial (in fact normalization constant) is given by Moreover, due to (70), we have Similarly to the case of a real-valued potential [51] (p. 95), one can see that is an entire function of ρ (see [51] (p. 95)). In this case, as a Birkhoff solution E(ρ, x), one can consider the Jost solution e(−ρ, x), Im(ρ) > 0, and hence from (72) we obtain Example 4. According to Remark 8, the normalization constant associated with the unique eigenvalue of the operator from Example 3 is and, in particular, for q 3 (x), we have Example 5. With the aid of Remark 8, an approximate value of the normalization constant for the eigenvalue λ 0 from Example 1 is obtained

Numerical Algorithm
The approximate solution of the direct problem can be performed with the following steps.

4.
To locate the eigenvalues, find the non-real poles of the function Ω(ρ, x), which is equivalent to finding zeroes of the function e(ρ) in the unit disk in terms of z. This can be achieved with the aid of the argument principle theorem. In particular, in the present work, we compute the change in the argument along rectangular contours γ.
If the change in the argument along γ is zero, consider another contour. Otherwise, subdivide the region within the contour until the desired accuracy is attained. Note that for a sufficiently large N, zeros of e N (ρ), approximate the square roots of the eigenvalues of the problem arbitrarily closely. The proof is analogous to that in [54] and is based on the Rouché theorem from complex analysis.

5.1
For simple poles, use Remark 8 to obtain the normalization constants.

5.2
Otherwise, for higher multiplicities, solve the linear system of Equation (65) for the coefficients b n k , n k = 0, 1, . . . , m k − 1 computing A j,n and D j defined in Equations (66) and (67) for several values of x j .

Inverse Problem
In order to reconstruct the potential in (1) from the scattering data, it is convenient to introduce the function [39] where η is a number satisfying the inequalities 0 < η < ε 0 (ε 0 s, defined in Section 3.2), and the function Remark 9. Hereinafter, we use the notation for 0 < η < ε 0 where L η represents a line parallel to the real axis crossing iη.
The kernel A(x, t) and the function f (x) satisfy the following Gel'fand-Levitan (G-L) equation [39] (Theorem 10.1) 4.1. Infinite Linear Algebraic System for Coefficients a n (x) Following [38], from the G-L Equation (78), we deduce the following system of linear algebraic equations for the coefficients a n (x) from the series representation (23).

Theorem 2.
The complex-valued functions a n (x) satisfy the equations Proof. Substitution of the series representation (23) into (78) leads to the equalities where the change in the order of summation and integration is justified by the general Parseval identity [47] (p. 16).
We have A(x, x + u), e −u/2 L n (u) (81) is equivalent to Multiplying the last equation by L m (s)e − s 2 and integrating this, we obtain Note that and ∞ 0 f (s + 2x + y)L n (y)e − y 2 dy = f n (2x + s).

Expressions for f m (x) and A mn (x)
It is convenient to regard the functions f m (x) and A mn (x) as a sum of the components corresponding to the continuous f m,c (x), A mn,c (x) and discrete spectra f m,d (x), A mn,d (x), and simplify these expressions with the aid of the formula ( [48], Formula 7.414 (6) The continuous and discrete components for the function f m (x) have the form and For the function A mn,c (x), we have and for A mn,d (x), we use (84) to obtain Remark 10. When an eigenvalue ρ 2 k is simple and the corresponding normalization polynomial X k (x) is just a normalization constant c k , expressions (86) and (88) can be written in the form We illustrate the calculation of the functions (85)-(88) with some examples.

Example 6.
Consider the scattering function obtained in Example 2: in the strip 0 ≤ Im(ρ) < 1, with no discrete spectrum and thus no normalization polynomials. Let us compute the function φ s (x) defined by (75), where the line L η lies in the strip 0 < Im(ρ) < 1.
Since the function s 2 (ρ) is analytic in the strip 0 < Im(ρ) < 1, the value of the integral is independent of the choice of 0 < η < 1. Using Jordan's lemma to calculate the integral in (75), we obtain Now, computing the functions f m (x) and A mn (x) from Formula (85) and (87) and using the residue theorem, we obtain Thus, in the case of the potential q 2 (x), the system of Equation (79) can be written explicitly.

Example 7.
Consider the scattering function s 3 (ρ) from Example 3. It has two poles in the upper half-plane: at i and i tanh(1 + i). Hence, using the residue theorem, we find that Again, the corresponding system of Equation (79) can be written explicitly.
Example 8. Consider s 1 (ρ) from Example 1. To compute φ s (x), we consider the singularities of s 1 (ρ) in the upper half-plane. From the set D(s 1 ) (see Example 1), we have that s 1 (ρ) has an infinite number of isolated singularities at the points ρ k = ik 2 with k ∈ N \ {0} and a singular number ρ 0 ; see (46). Using properties of the gamma function, we obtain where c 0 is the normalization constant obtained in Example 5. Therefore, for x > 0, we have Hence, the function f (x) has the form and we obtain the functions f m,c (x) and f m,d (x) in terms of z 0 = 1 2 +iρ 0 1 2 −iρ 0 (see (9)) as follows and f n,d (x) = (−1) n+1 c 0 e ixρ 0 (z 0 + 1)z n 0 .
Thus, from Equations (92) and (93) we obtain Likewise, applying the residue theorem, we have Thus, as in the previous two examples, the system of Equation (79) can be written explicitly.
The cancellation of terms when summing up (92) with (93) is not incidental and is generalized below in Remark 12.
To calculate the integrals in functions f m and A mn in the case when the scattering function is given explicitly, we implement Jordan's lemma and the residue theorem considering the asymptotics (44). However, often the function s(ρ) is not given in a closed form but as a table of data-then, the following techniques can be useful to compute the integrals. First, we recall a widely used technique for the quadrature of highly oscillatory integrals through approximations of the Fourier sine and cosine transform. This is illustrated below in Example 16. A second option is a transformation of integrals in f m and A mn into integrals over a finite interval providing a certain advantage for its numerical implementation. This is illustrated below in Example 21.

Remark 11.
We mainly discuss the calculation of the functions f m . The calculation of A mn is analogous.
provided the series on the right-hand side is convergent; see [55] for some 0 < η < 0 . Following the approach from [56] (p. 236), denote The integrals on the right hand side (the Fourier cosine and sine transforms) are approximated by the corresponding sums where h and N are chosen to be sufficiently small and large, respectively. 3.
Transform the line L η into a circle centered at −2η 1+2η of radius 1 2+η with the aid of the formulas This enables us to consider the integral in (85) in the form reducing the integration to a finite interval.

Stability of the System and Its Solution
Consider the truncated system (79): Denote its solution as U M = a M m M m=0 . In the following two theorems, we prove the unique solvability of (100), the convergence of its solution to the exact one as well as its stability. Proof. Since { f m (2x)} ∞ m=0 ∈ 2 and {A m,n (x)} ∞ m,n=0 ∈ 2 ⊗ 2 and we look for {a m (x)} ∞ m=0 ∈ 2 , the assertion of the theorem for the truncated system follows directly from the general theory presented in [57] (Chapter 14,  §3). Proof. Note that the truncated system (100) coincides with that obtained by applying the Bubnov-Galerkin procedure to the G-L Equation (78) with the orthonormal system of Laguerre polynomials in L 2 ([0, ∞); e −x ); see [58] ( §14). Let I M denote the (M + 1) × (M + 1) identity matrix, L M = {A m,n (x)} M m,n=0 be the coefficient matrix of the truncated system and R M = { f m (2x)} M m=0 the right hand side of (100). Following [58] ( §9), consider a system called inexact where Γ M is an (M + 1) × (M + 1) matrix representing errors in the coefficients A m,n , and δ M is the column vector representing errors in the coefficients f m . Let V M be a solution of the non-exact system. The solution of the Bubnov-Galerkin procedure is said to be stable if there exist constants c 1 , c 2 > 0, such that for Γ M ≤ r and arbitrary δ M the non-exact system is solvable, and the following inequality holds Now, since in the case under consideration, the inequality (102) is true (see [58] (Theorems 14.1 and 14.2)), the approximate solution is stable.

Algorithm to Recover the Potential
Given a scattering data set J as in Definition 2, the algorithm to recover q(x) consists of the following steps.

1.
Compute the functions f m (x) and A mn (x) with the aid of (85)-(88).

2.
Solve the truncated system of linear algebraic Equation (100) to obtain the coefficient a 0 (x).

Numerical Examples
We implemented the algorithms proposed in Sections 3.4 and 4.4 to solve the direct and inverse problems, respectively, with machine precision and with the aid of Matlab2021. Several examples are discussed, some of which have been introduced in previous sections.

Direct Problem
In this subsection, we discuss the computation of the scattering data, based on the series representation of the Jost solution (8). We deal with the approximate solution obtained by truncating the series (36).
The computation of the coefficients a n (x) is performed with the aid of the recurrent integration procedure from [27].
First of all, we discuss the choice of the number N in (36). Below, we show that a satisfactory accuracy is attained for a relatively small N (from several units to several dozens), and a reliable indicator can be used to choose an appropriate N.
In the case of simple singular numbers ρ k , the norming constants can be computed with the aid of (73): Another possibility consists of using (74) in the form where [m, n] e N (ρ) (−ρ k ) stands for the Padé approximant of e N (ρ) at ρ = −ρ k . This can be achieved when the accuracy of this rational approximation in the upper half-plane is satisfactory, i.e., when one has a suitable small value of max [m, n] e N (ρ) (ρ) − e N (ρ) in a sufficiently large region in the upper half-plane of the complex variable ρ.
A reliable algorithm to compute derivatives of (36) in (104) is proposed in [27].
To obtain the scattering function (42) in the strip |Im(ρ)| < ε 0 we consider two options depending on how the computation of the Jost function is performed for ρ in the lower half-plane. The first one uses provided [m, n] e N (ρ) (ρ) extends e N (ρ) analytically onto a certain strip in the lower half ρ-plane. A second option for the computation of s(ρ) is where the expression (36) is calculated at points ρ of a parallel line sufficiently close to the real axis and contained in the lower half ρ-plane. Remark 13. The notation for the approximate Jost solution (Jost function) may contain two indices, k and N: e k,N (ρ, x) (e k,N (ρ)), where k denotes the solution associated with the Schrödinger equation with the potential q k (x) and N is the parameter from (36).

Example 9.
Consider the potential q 2 (x) from Example 2. We present the indicator ε N in Table 1 for different values of N in (103).    Table 2 presents the maximum absolute and relative errors of the approximate Jost function e 2,N (ρ(z)) for z ∈ D for different values of N.   The distribution of the absolute and relative errors of the approximate Jost function is presented in Figure 3 and Figure 4 (respectively), where the maximum absolute error is 1.98 × 10 −14 and the maximum relative error is 3.17 × 10 −13 .  Furthermore, a good approximation of the derivative of the Jost function becomes essential for the argument principle algorithm performance. This is necessary to obtain the eigenvalues as the squares of non-real zeros of the approximate Jost function. In Figure 5, we illustrate de 2,30 (ρ(z)) dz , and Figures 6 and 7 depict the distribution of the absolute and relative errors, respectively. The maximum absolute error is 9.6 × 10 −13 and the maximum relative error is 7.21 × 10 −13 .  To find the singular numbers, we consider the circle {z ∈ C : |z| = 1} (real axis in ρ) and a cubic spline interpolation of the approximate Jost function (N = 30). For the spline interpolation, we use the Matlab routine csapi. To locate the zeros of the spline, we use slmsolve from the Shape Language Modeling (SLM) toolbox, version 1.14 by John D'Errico [59], available for Matlab2021a. The value ρ 1 = −1.000000000000003 was obtained with an absolute error of 3.11 × 10 −15 . Additionally, the argument principle algorithm applied to e 2,30 (ρ(z)) in D discarded any eigenvalue of the problem (non-real zero ρ).
In Table 3, we computed the maximum absolute and relative errors of the Padé approximant [1,1] Table 4, we confirm that this approximant satisfactorily extends the Jost function to a desirable strip in the lower half-plane (the strip is related to the one needed for the calculation of the scattering function s 2 (ρ)).  To obtain s 2 (ρ) numerically on the strip 0 < Im(ρ) < ε 0 = 1, we use the truncated series e 2,30 (ρ) and the Padé approximant [1,1]  as the most suitable option to avoid the appearance of Froissart doublets. Indeed, the use of the Padé approximants when there is no available information about the smoothness of the function to be approximated is challenging. Some publications propose modified algorithms [60], even using the Toeplitz matrix theory with many numerical implementations in Maple, Wolfram Mathematica (see [61]) or Matlab (see [62]). For the purposes of this paper, it is sufficient to use only the information obtained from the truncated series e N (ρ) and the argument principle algorithm to construct the approximant. Consider the number K of zeros counting multiplicities of the approximate Jost function e N (ρ) (singular numbers being calculated using the argument principle algorithm) located inside D as the degree of the polynomial in the numerator in the Padé approximant. Recalling that, in most cases, an accurate Padé's approximation is obtained on the diagonal approximant types for analytical functions, it is reasonable to choose the Padé approximant as [K, K] e N (ρ) .

Example 10.
Consider the potential q 3 (x) from Example 3. The approximate Jost function e 3,N (ρ) is computed in the strip 0 ≤ Im(ρ) < ε 2 = 1 for several values of N. In Table 5, the maximum absolute error of the approximate Jost function is presented.  Similarly to the previous example, a search for real singular numbers was performed; however, none were detected. Subsequently, the argument principle algorithm located a non-real singular number in D, with the value z 1 ≈ −0.386709149322063 − 0.105221869864471i (ρ 1 ≈ −0.271752585319512 + 1.083923327338694i). Its absolute error is 8 × 10 −15 . The contour refinement is not a concern, since the performed algorithm from [54] is based on the argument principle algorithm followed by several Newton iterations.
Additionally, the Jost function was extended to the strip |Im(ρ)| < ε 0 = 1 through the Padé approximant [1,1]  Next, an approximate value of the normalization constant corresponding to ρ 1 was computed with an absolute error of 2.8 × 10 −9 . Finally, we calculate the scattering function by The maximum absolute error of the approximation of s 3 (ρ) in R 1 is 1 × 10 −9 (see Figure 8). . Example 11. Consider the potential q 1 (x) from Example 1. Table 6 shows the parameter ε N for some values of N. Note that the approximation of the Jost function in this example requires more terms in the series representation than in previous examples. To control the accuracy of the approximation, in addition to the parameter ε N , one can use the asymptotic relation for the Jost function from [24] (p. 105), where ω(x) = − 1 2 ∞ x q(s)ds. This relation is valid for q(x) with first and second summable derivatives. Figure 9 depicts the Jost function computed with N = 98 and the singular number ρ 0 ≈ 1.784065846059995 + 0.608788673578742i. Figure 10 shows the fulfillment of the asymptotic relation (108), namely the graph of e 1,98 (ρ) − ω(0) iρ + q(0) (2iρ) 2 − ω 2 (0) (2iρ) 2 , which tends to 1 when |ρ| → ∞.
The eigenvalue is computed numerically as a zero of the exact Jost function with the aid of Wolfram Mathematica v.12 (Wolfram Research, Inc., Champaign, IL, USA) λ 0 ≈ λ * 0 := 2.8122672899483 + 2.1722381890043i. This "exact" eigenvalue is compared with the approximation 2.812267289948449 + 2.172238189004328i obtained as the square of the approximate ρ 0 . The absolute error is 1.52 × 10 −13 .
For the numerical calculation of the analytic extension of e 1,98 (ρ) onto the strip − 1 2 < Im(ρ) < 0, it is not possible to consider the Padé approximant [1, 1] e 1,98 (ρ) . This does not approximate e 1,98 (ρ) accurately even in the upper half-plane of the complex variable ρ. Using the Padé approximant [7,7]  To compute the scattering function s 1 (ρ) on a line parallel to the real axis contained in the strip |Im(ρ)| < ε 0 = |Im(ρ 0 )| ≈ 0.608788673578742i, Formula (106) was used. The function e 1 (ρ) is represented by (36) for ρ on a line in the lower half ρ-plane parallel and sufficiently close to the real axis. Having calculated these series representations for the functions involved in s 1 (ρ), we compute with a maximum absolute error 1.45 × 10 −7 along the line L η=0.1 (see (77)).
In this example, we obtain a satisfactory accuracy in the calculations of the scatterin data set using the expression (36) alone and the derivatives required by (104).
The Jost solution is not available in a closed form. In order to check the validity of the numerical calculation of the coefficients a n (x) for e 4,N (ρ.x), we consider the indicator ε N (Table 7).  Figure 11 depicts the Jost function computed with N = 137 and the approximation of the singular number ρ 1 ≈ 1.416695330664399 + 0.634534798062634i, with its square belonging to the box B. Additionally, Figure 12 shows the fulfillment of the asymptotic relation (108).  The normalization constant c 1 is calculated using (104), Finally, the scattering function is approximated by e 4,137 (−ρ) e 4,137 (ρ) . Now, take R = 30 in the potential q 4 (x). In this case, two boxes localizing the only two eigenvalues λ 1 and λ 2 were obtained in [44] 55 46 i. Table 8 provides the values of the indicator ε N for several values of N.  Note thatλ 1 ∈ B 1 andλ 2 ∈ B 2 for N = 200. Finally, the normalization constants are calculated using e 4,200 (ρ) in (104), Although, in this example, more powers for the series representation of the Jost function were used, the method proved to be applicable to obtaining the scattering data set without any additional informatio. The good accuracy achieved is confirmed by the ability t use the scattering data obtained as input data to solve the inverse scattering problem to recover the potential q 4 (x) with R = 30 below in Example 22.

Inverse Problem
In the present section, we discuss the accuracy, convergence and stability of the proposed method for solving the inverse scattering problem.

Remark 15.
By q k,M (x), we denote the approximation of the potential q k (x) (k = 1, 2, 3, 4, 5) obtained by solving the truncated system (100) with the sum up to M, i.e., with M + 1 equations.
We shall recover the potential q 3 (x) = −2 sech 2 (x − 1 − i). The system (100) of linear algebraic equations for this example is obtained in a closed form (see Example 7). For a different number of equations in the truncated system, we obtain a solution symbolically by using the Matlab routine solve. The potential q 3 (x) is recovered from (38). Figure 13 presents the recovered potential in each case. The corresponding absolute and relative errors are presented in Figure 14 and Figure 15, respectively. Note that a high accuracy is attained even in the case of a very reduced number of equations in the truncated system. Moreover, a very fast convergence of the method can be appreciated.

Example 14.
Consider the scattering data J = {s 2 (ρ)} from Example 2. As was shown above (Example 6), the system (100) for this example can be written explicitly. Again, when solving the corresponding truncated system for different values of M we observe a fast convergence and remarkable accuracy even for small values of M (see Table 11 and Figure 16).  Figure 16. Exact and computed potential q 2,6 (x).

Example 15.
Consider the closed form of the scattering function s 1 (ρ) from Example 1. We compute functions f m,c (x) and A mn,c (x) using the first option from Remark 11. Some poles and residues are given in Table 12 (computed with the aid of the package Numerical Calculus of Mathematica v.12). Note that the absolute value of the residues decreases considerably as the poles move away from the origin on the imaginary axis. This allows us to use a small number of poles for the calculation of the functions f m,c (x) and A mn,c (x).
The convergence of the method in this case results to be slower; see Figure 17, although a satisfactory accuracy is attained for M = 9.

Stability of the System
Since the stability of the method was proved in Theorem 4, we are able to work efficiently with noisy scattering data. First, we consider the natural noise arising from the numerical implementations of the last two procedures in Remark 11, i.e., calculation of the approximate matrix in (100) from the scattering function s(ρ) given in a closed form. Another situation considered in this subsection is the recovery of the potential from a uniformly noisy scattering function.

Remark 17.
In the last step of the algorithm from Section 4.4, for recovering q with the aid of (38), the coefficient a 0 needs to be differentiated twice. This was performed by interpolating a 0 (x) with a quintic spline through the Matlab routine spapi and a posterior differentiation with the Matlab command fnder.

Example 16.
Let us consider the scattering data from Example 2. The recovery of the potential q 2 (x) from the exact scattering function s 2 (ρ), obtained by using approximate functionsf m (x) andÃ mn (x) in the truncated system (100), is presented. The computation of functions f m (x) and A mn (x) requires numerical integration along the line L η=0.5 (see (77)). For this purpose, the last two procedures in Remark 11 were applied. Method 1. The second option in Remark 11 is implemented. With the scattering function (91) at points ρ = σ + 0.5i and σ = −(k + 1/2)h for k = 0, 1, . . . , N(x), where N(x) = 55000/x and h = 0.145454545, the calculation of the Fourier transforms in (98) is carried out. In Table 13, the maximum absolute error off m (x) is presented for 4 values of the parameter m. Now, we computeÃ mn,c (x) using the same numerical integration method with parameters N(x) = 5500/x and h = 0.127272727. Table 14 shows the maximum absolute error ofÃ mn (x) for parameters m, n = 0, 1, 2, 3.  The system (100) constructed withf m (x) andÃ mn (x) is solved numerically in Matlab for several values of M. Maximum absolute and relative errors of the approximation of the potential q 2,M are shown in Table 15.  Figure 18 presents the absolute value of the recovered q 2 potential from 4 equations in (100). Method 2. Now, we compute the approximate functionsf m (x) (see Table 16) andÃ mn (x) (see Table 17) following the third procedure in Remark 11.  In Table 18, the absolute error of the recovered potential for some values of M in (100) is presented. Both methods (procedures 2 and 3 from Remark 11) illustrated in the above example have proven to be suitable for calculating the functions f m and A mn from a table of values for the s 2 (ρ). Nevertheless, it is worth mentioning that although the first method (procedure 2) produced slightly more accurate results, this approach might be sensitive to the choice of the N(x) and h parameters, whereas the second method (procedure 3) only requires the implementation of trapz, the Matlab integration routine on a dense set of points defined in the interval (0, 2π). Hence, for the purposes of this paper, it is sufficient to consider procedure 3 from Remark 11 in the following examples, so as to obtain satisfactory approximations of f m and A mn .
As expected from the results of Example 14 for this potential, the numerical method for recovering the potential q 2 (x) converges very fast. Indeed, an acceptable approximation of q 2 (x) is achieved with only four equations in this case, where an inexact matrix in the linear system (100) is considered. In fact, the difference between the approximate and the exact potential presented in Figure 18 is indistinguishable.
In the following examples, a noisy scattering function with a uniformly distributed noise ε(ρ) added to the rand routine of Matlab is considered.

Example 17.
Consider the scattering function s 2 (ρ) and denote the noisy scattering function byŝ 2 (ρ) := s 2 (ρ) + ε(ρ). Here, ε(ρ) is ±5% uniformly distributed complex-valued noise (the percentage of the noise is applied pointwise to the modulus and argument of the value of s 2 (ρ)). The maximum absolute error ofŝ 2 on the line L η=0.5 is 2.46 × 10 −1 . The potential was recovered using five equations with a maximum absolute error of 5.2 × 10 −1 . The real and imaginary parts of the potential and the absolute error of its recovery are shown in Figure 19. Despite the noise thatŝ 2 (ρ) produces in the matrix of the system (100), the method recovers the shape of the potential q 2 with reasonable accuracy.
The maximum absolute error ofŝ 3 on the line L η=0.5 is 1.75. The potential was recovered using eight equations with a maximum absolute error of 8.6 × 10 −1 . The real and imaginary parts of the potential as well as the absolute error of its recovery are shown in Figure 20. Although, in this case, the absolute error ofŝ 3 (ρ) is larger, the shape of the recovered potential is still quite close to that of the exact one.

In-Out
In this subsection, we consider the results obtained in Section 5.1 as input data for the inverse problem.

Example 19.
We use the approximate scattering function s 3 (ρ) from Example 10 calculated by (107). Particularly, the form in which it is given allows for us to approximate functionsf m,c (x) andÃ mn,c (x) with the aid of the numerical calculus of residues, i.e., the first procedure in Remark 11 (see Tables 19 and 20).  The potential q 3 (x) was recovered with an absolute error of 1.8 × 10 −5 in the interval (0, 15) using 8 equations.
Example 20. Consider the approximate scattering function s 1 (ρ) from Example 11. The coefficient a 0 (x) was recovered using 14 equations with a maximum absolute error of 4.29 × 10 −3 , from which the potential was recovered with a maximum absolute error of 0.23, Figure 21. Example 21. Consider the approximate scattering data obtained in Example 9. The approximate functionsf m (x) (see Table 21) andÃ mn (x) (see Table 22) were obtained accurately enough to recover the potential (see Table 23).

Example 22.
Consider the potential q 4 (x) = 30i sin(x) exp(−x) introduced in Example 12. Using the results of the solution of the direct scattering problem from Example 12, we recover q 4 (x) using 20 equations with a maximum absolute error of 8.67 × 10 −1 (Figure 23). It is worth mentioning that the coefficient a 0 (x) is recovered with an absolute error of 2.28 × 10 −2 ( Figure 24). The error is calculated and compared with the solution of the Cauchy problem a 0 (x) − a 0 (x) = q 4 (x)(a 0 (x) + 1), for a sufficiently large value of b > 0, obtained using ode45 routine of Matlab2021a. This is a case where closed formulas for the scattering data set are unavailable. Therefore, the In-Out procedure confirms a satisfactory accuracy in the solution of both the direct and inverse scattering problems. Example 23. Consider the singular potential q 5 (x) = exp(−2.5x) x − π 2 1/3 .
In Table 24, we present the parameter ε N for different values of N in (103). Using data from Table 24, we computed the scattering data with N = 45. No eigenvalue was detected, so the scattering data set consists of the scattering function approximated by the expression s 5 (ρ) ≈ e 5,45 (−ρ) e 5,45 (ρ) , ρ ∈ R.
Using this scattering data set to solve the inverse problem, we obtained the coefficient a 0 (x) as shown in Figure 25. The maximum absolute error resulted in 1.9 × 10 −4 . The potential is recovered as shown in Figure 26. The corresponding absolute error is presented in Figure 27. Indeed, the maximum absolute error is 9.82 × 10 −2 .  This example shows the applicability of the proposed algorithms to both the solution of the direct and inverse scattering problems in the case of non-smooth potentials.

Conclusions
An approach to solving the direct and inverse scattering problems on the half-line for the one-dimensional Schrödinger equation with an exponentially decreasing complexvalued potential is developed. It is based on a series representation of the Jost solution from [25], which is shown in the present work to remain valid in a non-selfadjoining case.
When solving the direct problem, this representation is used to calculate the scattering data set through a simple and efficient procedure, which includes a proposed algorithm for computing normalization polynomials (which are part of the scattering data set) by solving a finite system of linear algebraic equations for its coefficients.
When solving the inverse problem, the use of the series representation combined with the Gel'fand-Levitan equation reduces the problem to a system of linear algebraic equations for the series coefficients, and the knowledge of the first coefficient is sufficient to recover the potential.
The numerical results illustrate the remarkable accuracy of the proposed algorithms in solving both the direct and inverse scattering problems.