We consider first the single linear second-order equation
subject to the general two-point boundary conditions
We assume the functions p(x), q(x), and r(x) to be continuous on [a, b] and require that the homogeneous problem
have only the trivial solution, z(x) ≡ 0. Then by the Alternative Theorem (see Problem 1.2.5) the boundary-value problem (2.1.1) has a solution which is unique. Of course we could have required (2.1.1) to have a unique solution, and then it would follow that (2.1.2) has only the trivial solution. We shall describe and analyze the initial-value or shooting method for computing accurate approximations to the solution of the linear boundary-value problem (2.1.1).
Two functions y(1)(x) and y(2)(x) are uniquely defined on [a, b] as solutions of the respective initial-value problems
and
Here c0and c1are any constants such that
The fact that these problems have unique solutions on [a, b] is an obvious consequence of Theorem 1.1.1. The function y(x) defined by
satisfies a0y(a) – a1y’(a) = α(a1c0 – a0c1) = α and so will be a solution of problem (2.1.1) if s is chosen such that
This equation is linear in s and has the single root
provided that [b0y(2)(b) +b1y(2)′(b)] ≠ 0. However, if this quantity does vanish, then z(x) ≡ y(2)(x) would be a nontrivial solution of the homogeneous problem (2.1.2). Since this has been excluded, we see that the above construction of a solution of (2.1.1) is valid.
The initial-value or shooting method consists in simply carrying out the above procedure numerically; that is, we compute approximations to y(1)(x), y(1)′(x), y(2)(x), and y(2)′(x) and use them in Equations (2.1.4). To solve the initial-value problems (2.1.3) and (2.1.4) numerically, we first write them as equivalent first-order systems, say
and
Now any of the numerical methods discussed in Section 1.3 can be employed. Of course, the same method should be used to solve both systems (2.1.5a–b), and for convenience we will require that this be done on the net
so that xJ = b.
Let the numerical solutions of (2.1.5) obtained on the net (2.1.6) be denoted, in obvious notation, by the net functions
We assume that the scheme employed has an order of accuracy at least p and is stable. Then if the solutions of problems (2.1.5a-b) are sufficiently smooth, we are assured that for 0 ≤ j ≤ J,
provided the initial data used in the calculations also satisfy these relations. At the point x0= a the exact data can be used so that
To approximate the solution, y(x), of the boundary-value problem (2.1.1) we now form, by analogy with Equations (2.1.4),
where
We may expect, in light of the estimates (2.1.7), that |Yj –y(xj)| = O(hp).This is indeed true but, as we shall also show, there may be practical difficulties in actually computing an accurate approximation to y(x) for a fixed small h > 0.
Let us define the errors in the numerical solutions by
Evaluate Equation (2.1.4a) at x = xj and subtract the result from Equation (2.1.9a) to get, with the above notation,
However, a small calculation, using Equations (2.1.4b) and (2.1.9b), reveals that
Recalling the estimates (2.1.7), we have
and so for sufficiently small h the denominator in Equation (2.1.10b) cannot vanish. From Equations (2.1.10) we now conclude that
It is no more difficult to demonstrate that
so we also have accurate approximations to the derivative of the solution, for sufficiently small h.
The above result can now be summarized as follows.
THEOREM 2.1.1. Let the boundary-value problem (2.1.1) have a unique solution, y(x). Let the numerical solutions of the initial-value problems (2.1.5a–b), using a stable method with order of accuracy p on the net (2.1.6), be denoted by 
, v = 1, 2. Then the net functions
with sh as given by Equation (2.1.9b) and h sufficiently small f satisfy
provided the initial-value problems (2.1.5a–-b) have sufficiently smooth solutions. 
A particularly important case of the general boundary-value problem (2.1.1) occurs when a0 = b0 = 1 and a1 = b1 = 0 in the boundary conditions (that is, y specified at both ends). The error expressions (2.1.10) yield in this special case
However, the initial condition imposed on yj(2) is now, from Equations (2.1.8) with a1 = 0, just y0(2) = 0. Thus we have e0 = eJ = 0 and our approximate solution must be very accurate near both endpoints. Of course if the ratio yj(2)/yJ(2) becomes large in the interior, a < x, < b, this accuracy may be greatly diminished. In particular, if the fixed number y(2)(b) is small, then yJ(2) = y(2)(b) + O(hp), which depends upon h, may be very close to zero for some particular mesh spacing.
Returning to the general problem, we see that a similar loss in accuracy or error magnification may occur whenever
becomes very large. Of course since this ratio involves computed quantities, a practical assessment of this difficulty is possible.
Another very important source of difficulty can be loss of significant digits in forming the expressions in Equations (2.1.9a). Thus if Yj(1) and shYJ(2) are nearly equal and of opposite sign for some range of j values, we get cancellation of the leading digits. Unfortunately this situation can easily arise. For instance suppose that the solution of the initial-value problem (2.1.3a) grows in magnitude as x → b and that the boundary condition at x = b has b1= 0 (that is, y(b) = β is specified). Then if 
we have, from Equation (2.1.9b), approximately
and so Equation (2.1.9a) becomes
Clearly the cancellation problem occurs here for xj near b. Note that the solution y(1)(x) need not grow very fast, and in fact for β = 0, corresponding to the commonly-occurring boundary condition y(b) = 0, the difficulty is always potentially present. If either of the solutions, y(1)(x) or y(2)(x),changes rapidly near x = b, then these errors will not propagate very far from the endpoint.
The loss of significant digits can also occur with b1 ≠ 0 but is always easily observed. However, when employing an automatic digital computer one must remember to test for this effect as there may be no sign of it in a list of the leading digits of {Yj}. If the loss of significance cannot be easily overcome by carrying additional figures, that is, by using double-precision arithmetic, then the parallel-shooting techniques of Section 2.4 can be employed.
There is no difficulty in extending the results of this section to the case in which a0 = 0 but a1 ≠ 0 in the boundary conditions (2.1.1b). We pose this as Problem 2.1.1.
The initial-value method can easily be applied to problems which are much more general than problem (2.1.1). To illustrate such applications we shall consider a two-point boundary-value problem for a system of n coupled linear second-order differential equations. As we shall see, such problems can be solved in terms of the solution of n + 1 initial-value problems for the system. More general linear (and nonlinear) systems are treated in Section 2.3.
We write the system as
where y(x) and r(x) are n-vectors and P(x) and Q(x) are nth-order matrices whose elements are continuous functions of x. The two-point boundary conditions are
Here A0, A1, B0, and B1are nth-order constant matrices, α and βare constant n-vectors and again we assume that the homogeneous problem, that is, α ≡ β ≡ r(x) ≡ 0 in (2.1.11), has only the trivial solution, y(x) ≡ 0. An alternative theorem holds for such systems and so we know that a unique solution exists.
We now introduce n + 1 vector-valued functions y(v)(x), v = 0, 1,…, n, as the solutions of the corresponding initial-value problems
The vector I(v) is the vth-unit vector, that is, the vth column of the nth-order identity matrix I. By Theorem 1.1.1 it follows that each initial-value problem in (2.1.12) has a unique solution on [a, b]. The solutions of problems (2.1.12bv) form an nth-order matrix which we denote by
Then with an arbitrary n-vector, s = (s1s2, …, sn)T, define the vector
Now if s is a root of
then the function y(x) = y(x; s) is the solution of the boundary-value problem (2.1.11).
Proceeding formally, with Equation (2.1.13) in Equation (2.1.14), we see that ϕ(s) is linear in s and the system (2.1.14) has the unique solution
provided the indicated inverse exists. However, if the matrix
is singular, then its columns must be linearly dependent. That is, since the columns are B0y(v)(b) +B1y(v)′(b), v= 1, 2, …, n, there must exist scalars av, v = 1, 2,…,n, not all vanishing, such that
satisfies B0z(b) +B1z’(b) = 0 and so is a solution of the homogeneous boundary-value problem corresponding to (2.1.11). It has been assumed that this problem has only the trivial solution, z(x) = 0, and since
it follows that av = 0, v = 1,2,…,n. From this contradiction we conclude that Equation (2.1.14) indeed has the unique solution shown in Equation (2.1.15).
Numerical solutions of (2.1.11) are obtained by solving each of the n + 1 initial-value problems in (2.1.12) and forming, at the net points, the quantities corresponding to (2.1.13) and (2.1.15). The linear system (2.1.14) which must be solved is usually of relatively small order and so is not, in principle, a difficult numerical problem.
To estimate the accuracy of the initial-value method applied to the problem (2.1.11) we proceed as in the previous case (where n = 1). Thus we need not write the obvious first-order systems which are solved in place of (2.1.12). Rather we denote by 
, and s(h) the numerical approximations to y(xj), y(v)(xj), y(v)’(xj), and s, respectively. If a stable method with order of accuracy at least p is employed to solve the initial value problems then we may assume that
Now subtracting Equation (2.1.13) with x = xj from its numerical counterpart we get
Similarly from Equation (2.1.14) and its numerical equivalent we find that
where the matrix H is
Here E is the matrix with columns 
, v = 1, 2, …, n, and 
has the columns 
, v = 1, 2,…,n. Thus for sufficiently small h the matrix H is nonsingular, since the matrix (2.1.16) is nonsingular (see Problem 2.1.3), and we conclude that
just as in the case with n = 1. A similar estimate applies for the first derivatives.
Combining these observations we may easily construct a proof of the following.
THEOREM 2.1.2. Let the boundary-value problem (2.1.11) have a unique solution, y(x). Let the numerical solutions of problems (2.1.12a, bv), using a stable method with order of accuracy p on the net (2.1.6), be denoted by 
, 
, v = 0, 1,…, n. Define the net functions
For sufficiently small h a unique vector 
is defined as the solution of the linear system
With this root the net functions {Yj(s(h))} and {Vj(s(h))} satisfy
provided that the initial-value problems (2.1.12) have sufficiently smooth solutions.
Although this theorem is similar to Theorem 2.1.1 and the error estimates, valid as h → 0, are gratifying, it should not be inferred that the practical difficulties are as simple as they were in the case n = 1. Combining Equations (2.1.17) and (2.1.18), we see that if H is ill-conditioned for finite values of h we may lose accuracy. The cancellation difficulty in forming (2.1.13) is now even more important. It can occur, for instance, if any fixed corresponding components of 
and 
for fixed v1 ≠ v2are nearly equal, of opposite sign, and in magnitude larger than the other corresponding components for v ≠ v1, v2. But again some of the sources of error may be eliminated by using the parallel-shooting methods of Section 2.4.
2.1.1 Replace the condition detA0≠ 0 in Equations (2.1.11b) by the condition rank [A0, A1] = n and show how to solve the resulting boundary-value problem (2.1.11).
2.1.2 Prove the Alternative Theorem for the linear second-order system (2.1.11). [HINT: Use the reduction to an initial-value problem for one system which has a unique solution.]
2.1.3* Show that A + B is nonsingular if A is nonsingular and the elements of B are sufficiently small [see Isaacson and Keller (1966), Theorem 1.5, Chapter 1].
2.1.4 Use the above result to show that Equation (2.1.19) has a unique solution s(h) for sufficiently small h. [HINT: Relate the coefficient matrix to H in Equation (2.1.18b).]
2.1.5 Use the initial-value method to solve the boundary-value problem
If the linearly-independent solutions
are employed, show that, for small ||,
Discuss the possible computational difficulties caused by small ||.
2.1.6 Use the initial-value method to solve the boundary-value problem:
If the linearly-independent solutions
are employed, show that
Explain the difficulty in carrying out the numerical solution for small ||.
We now consider second-order equations, which may be nonlinear, of the form
subject to the general two-point boundary conditions
The function f(x, y, z) in Equation (2.2.1a) will be assumed to satisfy the hypothesis of Theorem 1.2.2. This assures us that the boundary-value problem (2.2.1) has a solution which is unique.
The related initial-value problem that we consider is
where c0 and c1 are any constants such that
The solution of this problem, which we denote by u = u(x; s), will be a solution of the boundary-value problem (2.2.1) if and only if s is a root of
Under the previously-stated conditions on f, it is shown in the proof of Theorem 1.2.2 that the function ϕ(s) has a positive derivative which is bounded away from zero for all s, and so Equation (2.2.3) has a unique root for any value of β.
The initial-value or shooting methods we now study essentially consist of iterative schemes for approximating the root of Equation (2.2.3). In particular we shall employ schemes of the type indicated in Section 1.4. These all require evaluations of the function ϕ(s) for a sequence of values of s. This is done by means of numerical solutions of a sequence of initial-value problems of the form (2.2.2). [Of course if f(x, y, z) is linear in y and z, then ϕ(s) is linear in s and the root of ϕ(s) = 0 could be determined by only two evaluations of ϕ(s). This is what was done in Section 2.1. However, a linear equation can also be solved by iteration and so the present procedures also apply to linear boundary-value problems.] Obviously we should try to employ very rapidly converging iteration schemes as this would reduce the number of initial-value problems that must be solved. We shall discuss the use of several different schemes, one of which is Newton’s method and is particularly well-suited for the current class of problems. First we must consider some effects of using numerical approximations to u(b; s) and u’(b; s) in attempting to evaluate ϕ(s).
To solve the initial-value problem (2.2.2) numerically, we write it as a first-order system, with v ≡ u′;
If we write the solution of this system as
function ϕ(s) as defined in (2.2.3) becomes
Again, any of the methods in Section 1.3 can be used to solve the initial-value problem (2.2.4) numerically. We use the net in (2.1.6) and denote the numerical solution, for each value of s, by
If the numerical method has order of accuracy p and is stable, and if the solutions of (2.2.4) are sufficiently smooth, then
In terms of the numerical solution we define the function
and from Equations (2.2.6) and (2.2.5) it follows that
Thus in general we cannot hope to evaluate ϕ(s) exactly but rather we get an O(hp) approximation to it.
It is clear, from the above error bounds, that the solution of the boundary-value problem (2.2.1) can be approximated to an accuracy at least O(hp) if the root s = s* of Equation (2.2.3) is known. However we can at best hope to approximate s* to an accuracy of O(hp) by virtue of the estimate (2.2.8). But then the solution of the boundary-value problem will also be approximated to this same accuracy by means of the following.
LEMMA 2.2.1. Let u (x; s) and v(x; s) denote the solution of the initial-value problem (2.2.4) with f ≡ (v,f(x, u, v))T satisfying the hypothesis of Theorem 1.1.1. Let {Uj(s)}, {Vj(s)} denote the numerical solution of (2.2.4) on the net in (2.1.6) using a stable, pth-order-accurate scheme. Then if the solutions of (2.2.4) are sufficiently smooth, we have for any values s and t
Proof. From the hypothesis it follows by Theorem 1.1.1 that u and v are Lipschitz-continuous in s. But then
where we have used the estimate (2.2.6b), which is valid, and K’ is the Lip-schitz constant for u. A similar result holds for v.
If in (2.2.9) we set s = s* and |s*– t | = O(hp) we see that both the solution of the boundary-value problem (2.2.1) and its first derivative can be computed to within O(hp) if s* can be approximated to this accuracy.
To approximate the root of Equation (2.2.3) we write the equation in the form
and consider the functional iteration scheme: s0= arbitrary,
The fact that the sequence {sv} converges to the root under quite general circumstances is the content of the following.
THEOREM 2.2.1. Let f(x, y, y’) in Equation (2.2.1a) satisfy the hypothesis of Theorem 1.2.2, and in addition for some positive constant N,
Then the function ϕ(s) in Equation (2.2.3) has a continuous derivative, 
, which satisfies
where
For any s0 and fixed m in
the iterates in Equations (2.2.10b) converge to the root s* of Equation (2.2.3). In particular, with the choice m = 2/(Г + γ), the iterates satisfy
Proof. We first note that Theorem 1.1.2 is applicable to the system (2.2.4). Thus the functions ξ(x) ≡ ∂u(x;s)/∂s and ξ’(x) ≡ ∂u’(x;s)/∂s, where u(x; s) is the solution of (2.2.2), are continuous functions of s for all x ∈ [a, b]. But then
is also continuous in s.
As in the proof of Theorem 1.2.2 we have the variational problem for ξ(x),
where now
It has already been shown [see the inequalities (1.2.5)] that since a0 ≥ 0 and a1 ≥ 0,
and so clearly, with the notation in (2.2.11), ≥ γ. To obtain an upper bound we observe that now
It easily follows (see Problem 2.2.1) that ξ(x) ≤ η(x) and ξ’(x) ≤ η’(x), where η(x) is the solution of
But this linear problem with constant coefficients has the solution
and so ≤ Г follows, establishing (2.2.11a).
We now use the mean-value theorem to deduce that
Thus g(s) in Equations (2.2.10) satisfies a Lipschitz condition with constant λ given by
For any m in 0 < m < 2/Г, it follows that λ < 1, and thus Theorem 1.4.1 is applicable (with ρ = ∞) to the scheme (2.2.10b). With m = 2/(Г+ γ) we obtain
and part (d) of Theorem 1.4.1 yields (2.2.12).
One of the nice features of this result is that it determines a specific value for m, in terms of known data and the two bounds M and N, which insures convergence. However, if we consider for the moment the case with a1 = b1 = 0, a0 = b0 = 1 (that is, function specified at both endpoints), then Equation (2.2.11b) yields
From this we see that the convergence factor, (Г –γ)/(Г +γ), in (2.2.12) may be very close to unity if the interval (b – a) is “long” or the constants M or N are “large.” In such cases many iterations may be required to obtain an accurate approximation to the root. However, we shall consider more rapidly convergent schemes shortly.
Let us now consider the effects on the above procedure which result from the fact that ϕ(s) cannot be evaluated exactly. That is, we employ the iteration scheme in Equations (2.2.10), using the numerical solutions (2.2.6) of (2.2.2) to evaluate the function ϕ(s). Recalling Equation (2.2.7), the actual iteration scheme is then: s0 = arbitrary,
The basic result can be stated as follows.
THEOREM 2.2.2. Let f(x, y, y’) in Equation (2.2.1a) satisfy the hypothesis of Theorem 2.2.1. Let numerical solutions (2.2.6) of the initial-value problem (2.2.4) be computed as in Lemma 2.2A. Let the sequence {sv} be defined by Equations (2.2.13), with Φ(s) as defined in Equation (2.2.7). Then for any m in 0 < m < 2/Г we have, with y(x), the solution of the boundary-value problem (2.2.1), for some positive λ = λ(m) < 1,
Proof. We may apply Lemma 2.2.1 with t = sv and s = s*, the root of Equation (2.2.3), in which case u(x; s*) = y(x). The proof is then reduced to estimating |s* –sv|. But from the estimate (2.2.8), which is valid by hypothesis, we can write Equations (2.2.13) as
It has been shown in Theorem 2.2.1 that g(s) ≡ s – mϕ(s) is Lipschitz-continuous with constant λ < 1for m in 0 < m ≤ 2/Г. Now Theorem 1.4.2 can be applied, with 
, to deduce that
It should be noted that the convergence factor λ, in the above theorem, is independent of the net spacing h. Thus for each h one can find an integer v = v(h, p, λ) such that
and so, in principle, the total error can be made O(hp).
The virtue of rapidly converging iteration schemes for solving Equation (2.2.3) is quite clear and we recall, from Theorem 1.4.3, that Newton’s method may have arbitrarily small convergence factors, λ. Thus we indicate the proper way in which to employ this method in the present problem. The iterations {sv} are now defined by: s0 = arbitrary,
We have previously indicated how ϕ(s) in Equation (2.2.5) is to be evaluated, or rather approximated, by Φ(s), in terms of a numerical solution of the initial-value problem (2.2.4). We shall evaluate by means of a numerical solution of the variational problem corresponding to (2.2.4). This variational problem, or one equivalent to it, has already been employed in Theorem 2.2.1 to study the properties of . If we introduce
where u and v = u’ are the solution of (2.2.4), then formal differentiation in this initial-value problem yields
where
Thus we have, from Equation (2.2.5),
The numerical solution of the initial-value problem (2.2.15) should be evaluated on the net (2.1.6), along with the numerical solution of the initial-value problem (2.2.4) and by the same method. Then in many cases the coefficients p(x) and q(x) can be evaluated with little additional computation over that required to evaluate f(x, u, v). Thus one iteration in Newton’s method (2.2.14) requires the solution of two initial-value problems, but the amount of computation may frequently be less than that required for two iterations in the scheme (2.2.10).
Finally, we point out the obvious but frequently forgotten fact that “shooting” can be applied in either direction. This applies to linear problems as well. If the solutions of the initial-value problem grow from x = a to x = b, then it is likely that the shooting method will be most effective in reverse; that is, using x = b as the initial point. We call this procedure reverse shooting.
2.2.1 Let ξ(x) and η(x) satisfy on [a, b]:
If the constants M and N satisfy M ≥ 0and N ≥ 0, show that ξ(x) ≤ η(x) and ξ’(x) ≤ η’(x) on [a, b]. [HINT: “Solve” the differential inequality satisfied by ζ(x) ≡ η(x) – ξ(x). First derive the result that [ζ’ + Bζ] ≥ 0, where 
. Then it follows that ζ(x) ≥ 0 on [a, b] and so [ζ” + M ζ’] ≥ 0 which implies ζ’(x) ≥ 0.]
2.2.2* If f(x, u, v) is as in Theorem 2.2.1 and in addition ∂f(x, u, v)/∂u and ∂f(x, u, v)/∂v are Lipschitz-continuous functions of u and v, show that as defined in Equation (2.2.16) is a Lipschitz-continuous function of s.
2.2.3 Use Theorem 2.2.1, Problem 2.2.2, and Problem 1.4.7 to show that Newton’s method in (2.2.14) converges if |b – a| is sufficiently small and some later iterate in (2.2.10) is used as initial guess.
2.2.4 Show that the choice m = 2/(Г+ γ) in Theorem 2.2.1 is in general the best possible in light of the inequalities (2.2.11a).
by shooting. Try using several different uniformly-spaced nets (that is, different values of h), contracting maps, and Newton’s method.
2.2.6 Verify the applicability of the theoretical results to the problem formulated in Problem 2.2.5.
The application of initial-value techniques to rather general nonlinear boundary-value problems is immediate. First of all, the particular equation or set of equations can be replaced by an equivalent system of say n first-order equations of the form
The components, yj(x), of the n-vector yare either the original dependent variables or some derivatives of them or linear combinations of these functions. Assuming the original two-point boundary conditions to be linear, we can now write them in the form
Here A and B are constant square matrices of order n and α is a specified n-vector. In the indicated reduction to (2.3.1) we have assumed that the highest-order derivatives of each dependent variable which occurred in the original system of differential equations did not occur in the boundary conditions. The n boundary conditions in Equation (2.3.1b) will be independent if the matrix (A, B), which is n × 2n, has rank n. It should be noted that periodic boundary conditions are included as a special case in (2.3.1b). For purposes of our general discussion we associate with the boundary-value problem (2.3.1) the initial-value problem
Under appropriate smoothness conditions on f(x,u) we are assured, as in Theorems 1.1.1 and 1.1.2, that a unique solution of (2.3.2) exists on a ≤ x ≤ b, say
and it depends Lipschitz-continuously or even differentiably on the components of s. We now seek ssuch that u(s; x) is a solution of the boundary-value problem (2.3.1). This occurs if and only if sis a root of the system of n equations
(see Theorem 1.2.4). In particular cases the initial values in Equation (2.3.2b) may be chosen quite conveniently, for instance so that say p of the conditions in (2.3.1b) are automatically satisfied. Then only q = n – p undetermined components sj would occur in the chosen replacement of (2.3.2b). These components must then be determined in order to satisfy the q conditions in (2.3.3) that are not identically satisfied. This procedure can be applied whenever p independent conditions in Equation (2.3.1b) involve only p components of y(a) [see Problem 1.2.12].
Under the hypothesis of Theorem 1.2.6, which we assume to hold, we are assured that Equation (2.3.3) has a unique root, or equivalently that the boundary-value problem (2.3.1) has a unique solution. Furthermore the root of Equation (2.3.3) is also the fixed point of a contracting map which we write as
Here, since (A + B) is assumed nonsingular, we take
(We recall that the most general boundary conditions (2.3.1b) do not necessarily satisfy this condition and the computations in such cases may converge with a different choice for Q.) Our formal procedure is now clear. To show that the solution of the boundary-value problem (2.3.1) can be computed to any desired accuracy, we approximate the root of Equation (2.3.4) by iterations and use a sufficiently accurate such approximation to “solve” the initial-value problem (2.3.2) numerically. Of course, in carrying out the implied iterations, a sequence of numerical solutions of (2.3.2) is required. We indicate below more details in this convergence proof which is quite similar to that in Section 2.2, and then we discuss the practical aspects of actually computing an accurate approximation. Newton’s method is again advocated and we show that it is applicable in the present case.
Let one of the numerical methods of Section 1.3 be used to obtain a numerical solution of the initial-value problem (2.3.2) on the net (2.1.6). We denote the numerical solution at each net point xj, for each value of s, by the n-vectors
If the numerical method is stable and has order of accuracy p, and the solutions of the initial-value problem (2.3.2) are, as we assume, sufficiently smooth, then
Now, just as in Lemma 2.2.1, it follows that for any two sets of initial data sand t
If s is a root of Equation (2.3.3) all we need do, by the estimate (2.3.6), is to find a t such that |s– t| = 
, and then Uj(t) is accurate of order p.
From the proof of Theorem 1.2.6 we know that g(s) in (2.3.4) satisfies
Thus, as previously stated, g(s) is contracting [by Problem 1.2.11]. However, in terms of the numerical solution, we approximate g(s) by
and the iterations employed to obtain the root of Equation (2.3.4) are s(0) = arbitrary,
However, by (2.3.5), (2.3.4), and (2.3.8a), we have
Thus the iterates (2.3.8b) actually satisfy
But since g(s) is contracting, it follows that Theorem 1.4.2 is applicable. Then if s* is the root of Equation (2.3.4), the iterates in (2.3.8) must satisfy
Using this result with t = s(v) and s = s* in (2.3.6), we obtain what may be summarized as follows.
THEOREM 2.3.1. Let f(x, y), A and B satisfy the hypothesis of Theorem 1.2.6. Further, let f(x, y) be so smooth that the numerical solutions Uj(s) of the initial-value problem (2.3.2) computed by a stable pth-order method on the net (2.1.6) satisfy (2.3.5). Then the sequence s(v)defined in (2.3.8) is such that
where y(x) = u(s* ; x) is the unique solution of the boundary-value problem (2.3.1) and λ < 1 satisfies the inequality (1.2.18e).
Of course the iteration scheme (2.3.8) is frequently not of practical value (when, for example, 1 – λ is too small), and we would generally prefer a higher-order method for greater efficiency. For this purpose we recall Newton’s method for approximating a root of Equation (2.3.3) as follows: s(0) = arbitrary,
where ∆s(v) is the solution of the nth-order linear algebraic system
These iterates are well defined when the coefficient matrix, or Jacobian matrix, ∂ϕ/∂sis nonsingular for each s(v). If this can be established, then the question of convergence is relevant. However, in the present case, Theorem 1.4.3 becomes applicable when the Jacobian is nonsingular. Thus to justify Newton’s method applied to Equation (2.3.3) we need only show that ∂ϕ(s)/∂s is nonsingular.
Differentiating in Equation (2.3.3), we obtain the representation
Here the nth-order matrix W(s; x) = ∂u(s; x)/∂s is the solution of the variational system, obtained by differentiating in (2.3.2),
We may write Equation (2.3.10), recalling the assumption that (A + B) is nonsingular, as
But it is shown in the proof of Theorem 1.2.6 that
and so we can conclude that the Jacobian matrix is nonsingular (see Isaacson and Keller (1966), p. 16) for all s. Of course, the matrix (2.3.10) may very well be nonsingular even though (A + B) is singular. Thus, Newton’s method is frequently applicable when our present proof of the fact is not.
The calculational procedure for employing Newton’s method consists in solving n + 1 initial-value problems for nth-order systems at each iteration. The single nonlinear system (2.3.2) is integrated numerically with s= s(v), say, and at the same time and on the same net the n linear systems in (2.3.11) are integrated using the currently-computed values of Uj(s(v)) in place of u(s(v); xj) to evaluate the elements in ∂f(x,u)/∂u. Then s(v+1) is computed from (2.3.9) where, however, ϕ(s) in Equation (2.3.3) and ∂ϕ(s)/∂s in Equation (2.3.10) are approximated using the appropriate numerical solutions for u and W at xJ= b.
There may be great practical merit in applying Newton’s method to linear systems of differential equations. Sometimes the iteration schemes are less sensitive to cancellation errors that can disturb and even invalidate direct applications of shooting, as pointed out in Section 2.1. More important, however, is the fact that the matrix W(x), the solution of (2.3.11), is independent of sfor linear problems. That is, if f(x, u) has the form
where K(x) is an nth-order matrix, then clearly
and the variational problem (2.3.11) reduces to
Recalling Equation (2.3.10), we then see that ∂ϕ(s)/∂s is independent of s and hence the system of Equations (2.3.3) for the determination of s is a linear system. (This was shown before more directly, in Section 1.2.1.) Since W(x) ≡ ∂u(s; x)/∂s, and we have just shown that ϕ(s) is linear, so that ϕ(s) = [∂ϕ(s)/∂s]s + ϕ(0), it follows that Equation (2.3.3) can be written as
Thus by solving the n systems (2.3.13) the coefficient matrix is determined, and by solving the single system (2.3.2) with s = 0 and f as given by Equation (2.3.12), the inhomogeneous data are determined.
Now to solve a sequence of linear boundary-value problems, given by (2.3.1) and (2.3.12), in which the boundary conditions change, we need only use the new data A, B, and α in Equation (2.3.14), solve the resulting linear algebraic system for s and then solve the initial-value problem (2.3.2) using this value of s. If in addition the inhomogeneous term in the differential equation is to be changed [this is the term f0(x) in Equation (2.3.12)], we need one additional integration of (2.3.2) with s = 0, as above, to compute u(0; b). In many linear problems of interest it is only the inhomogeneous data that are to be altered and in such cases the above-indicated procedures are very efficient. It may frequently be worthwhile, in such cases, to compute W(x) with extra accuracy (that is, multiple precision).
Of course in most nonlinear problems of practical interest the hypothesis of Theorem 1.2.6 will not be satisfied. We have already seen, for instance, that the condition, (A + B) nonsingular, is rather special. Nevertheless, the initial-value method we have indicated is very frequently applicable though the proofs are not. The practical problem is to determine some matrix Q, now ≠ (A +B)–1, so that (2.3.4) is contracting; see for example Problems 2.3.4 and 2.3.5. Much more general results are indicated in Section 5.5 and, for instance, in Conti (1958) or Losota and Opial (1966). Of course if the initial unknown vector scan be reduced to a q = n – p-dimensional vector, as previously discussed, this should be done. Then the system of Equations (2.3.3) reduces to q equations, the corresponding Jacobian matrix in Equation (2.3.10) is of order q and the matrix W is n × q. Now only q + 1 systems of order n need be integrated for each iteration [see Problem 1.2.12]. Needless to say this reduction should be done for linear problems when applicable.
It may be a formidable task to apply any iteration scheme successfully when one does not have some reasonable estimate of the location of a zero of ϕ(s). Unfortunately there are no universal procedures for obtaining such estimates. But patience in searching for an acceptable first estimate is frequently rewarded. It is not unusual in practice to devote ninety percent or more of the computing effort to locating a neighborhood of a root. In many problems of scientific interest, as opposed to textbook problems, there are parameters in the equations or boundary conditions, and solutions are required over an entire domain of these parameters. But then when a solution has been obtained for one set of parameter values we can frequently make small changes in the values of the parameters and obtain neighboring solutions with ease.
This device of altering a parameter can be so effective that several techniques for solving nonlinear problems are based on introducing such a parameter. For instance we may introduce a single parameter t in such a way that for t = 1 the problem is the one whose solution is desired and for t = 0 it is some problem that is easily solved. Now if the solution depends continuously on t we may be able to proceed in small steps ∆t, starting from t = 0 as indicated above, to compute finally the solution for f=1. The study of these methods is briefly considered in Section 5.5.1. Let it suffice to say here that two very powerful techniques for proving existence theorems—the continuity and Poincaré continuation methods—are based on this idea [see Bers, John, and Schechter (1964), pp. 238–240, 285, Lasota and Opial (1966) and Section 5.5].
2.3.1 Show that the necessary and sufficient condition for the n linear boundary conditions in (2.3.1b) to be linearly independent is that rank (A, B) = n [here A and B are nth-order matrices while (A, B) is a matrix with n rows and 2n columns].
2.3.2 Show that a necessary and sufficient condition for rank (A, B) = n is that there be two nth-order matrices, say P1 and P2, such that the nth-order matrix AP1 + BP2 is nonsingular.
2.3.3 Use the results stated in Problems 2.3.1 and 2.3.2 to deduce that a necessary condition for the applicability of Newton’s method in solving (2.3.1), that is, in finding a root of Equation (2.3.3), is that the boundary Conditions (2.3.1b) be linearly independent.
2.3.4* If ∂f(x,u)/∂u is Lipschitz-continuous in u, the matrix
in Equation (2.3.10) is nonsingular and |ϕ(s(0))| is sufficiently small, then Equation (2.3.3) has a root if |b –a|is sufficiently small. Use the result in Problem 1.4.5 of Chapter 1 to indicate a proof of the above and estimate the magnitudes of |ϕ(s(0))|and |b –- a| in terms of the magnitudes of |f |,
∂f/∂u and the Lipschitz constant.
2.3.5 Formulate the boundary-value problem
in terms of a first-order system. Show that the conditions of Theorems 1.2.5 and 1.2.6 are satisfied if L < log 2. For L > log 2, show that a contraction is obtained by using, in place of Q = (A + B)–1 = I, the matrix
with any q satisfying
What happens when cos L = 0? Try actual computations in the above with varying values of L.
In the previous sections we have discussed some of the difficulties which arise in the shooting method. Of those discussed the best known, perhaps because it is the most common, is that caused by the growth of solutions of the initial-value problems which must be solved. We recall that this causes a loss in accuracy in attempts to solve the corresponding algebraic or transcendental system. [Technically we may say that the problem of finding a root of ϕ(s) = 0 is “ill-conditioned” or “not properly posed” in the sense of Section 3, Chapter 1, in Isaacson and Keller (1966).] As has been pointed out, greater accuracy in the computations may sometimes overcome this difficulty. This is not always practicable or desirable and so we shall describe variations of the shooting method that can frequently eliminate this difficulty. The finite-difference or integral-equation methods of Chapters 3 and 4 do not suffer from this particular growth problem, in principle.
There is another more striking phenomenon that can hamper the initial-value method when it is applied to nonlinear problems. We have already pointed out that many, if not most, boundary-value problems of interest do not satisfy the nice conditions imposed in results like our Theorems 2.2.2 and 2.3.1 which insure the applicability of shooting. In particular the functions f(x, y, y’) or f(x, y) may be unbounded in the x, y-domain of interest. Then solutions of the corresponding initial-value problems may become singular in [a, b] for some values of the initial data. If this occurs we cannot even carry out the integration from a to b in order to improve the initial data. Furthermore, the location of the singularity is frequently a very sensitive function of the initial data and thus great accuracy in the initial estimate of s(0)ors(0) may be order to start an iterative scheme ; see Problem 2.4.1. When this difficulty with singularities occurs it is quite apparent in the computations. The parallel shooting schemes presented below can frequently circumvent the trouble, or else a finite difference procedure can be employed. As we shall see, parallel shooting is actually a combination of difference methods and initial-value methods.
The variant known as multiple or parallel shooting is designed to reduce the growth of the solutions of the initial-value problems which must be solved. This is done by dividing the interval [a, b] into a number of subinter-vals, integrating appropriate initial-value problems over each subinterval, and then simultaneously adjusting all of the “initial” data in order to satisfy the boundary conditions and appropriate continuity conditions at the subdivision points. The variety of procedures afforded by this technique is limitless, and we consider only two rather obvious forms of it here. Of course the motivation in selecting subintervals and directions of integration is, as indicated above, to reduce the magnification of errors caused by growing solutions.
To illustrate we shall consider the boundary-value problem, on an interval a ≤ x ≤ b, for a, system of n first-order equations:
We divide this interval into J subintervals with the points xj, say
We refer to the interval xj – 1 ≤ x ≤ xj, as the jth subinterval and denote its length by
Now on each subinterval we introduce the new independent variable t, the new dependent variable yj(t) and the new n-vector functions fj(t, z) by
With these changes of variables the system (2.4.1a) becomes, on the jth subinterval,
The boundary conditions (2.4.1b) are now
Continuity of the solution of (2.4.1) at each interior point xj is expressed by the conditions
The J systems of n first-order equations (2.4.4a) can be written in the vector form
and the J sets of conditions (2.4.4b, c) can be written in the matrix-vector form
Here we have introduced the nJ-dimensional vectors
and the square matrices of order nJ
In (2.4.7) the identity matrix I is of order n and the blocks of zeros are also of order n.
Under modest smoothness conditions on f(x, y), the boundary-value problem (2.4.1) for n first-order equations on a ≤ x ≤ b can be shown to be equivalent to the boundary-value problem (2.4.5) for nJ first-order equations on 0 ≤ t ≤ 1. Clearly one parallel shooting procedure for solving (2.4.1) consists in applying ordinary shooting, as in Section 2.3, to (2.4.5). For this purpose we must solve the initial-value problem
for a first-order system of nJ equations. We try to pick the nJ-dimensional initial vector, s, such that
Of course Newton’s method is advocated for solving this system. To employ it we compute, along with the solutions of the initial-value problem (2.4.8a), the nJ-order matrix solution, W(s, t) ≡ ∂U/∂s, of the variational system
The iterates are then given by s(v + 1) = s(v) +∆s(v), where
The computations involved in carrying out the above procedure are not nearly as complicated as they might seem. This is due to the special form of F(t, Y) defined in (2.4.6). In fact, if we introduce the n-vectors sj such that 
, then the system (2.4.8a) is in fact
This is J first-order systems of n equations and each system is independent of the others. Thus these J initial-value problems can be solved independently of each other. In fact the solution of the set of initial-value problems (2.4.8) is ideally suited for computation on parallel computers (which are at present being designed). The same is of course true of the variational system (2.4.10) for the matrix W of order nJ. If we define the J matrices, Wj, of order n by
then W(s, t) = diag {Wj(sj;t)}. This follows since Wj = ∂uj/∂sj, and ∂ui/∂sj ≡ 0 if i ≠ j. The only coupling of the systems (2.4.8b) and (2.4.10b) for different values of j occurs in the algebraic problem of solving the linear system (2.4.11). We examine below some questions relating to the nonsin-gularity of the coefficient matrix in Equation (2.4.11). In particular the explicit form of the inverse of this matrix can be deduced under special circumstances (see Problem 2.4.3).
We first observe that the nJ-order matrix (P +Q) is nonsingular if the nth-order matrix (A + B) is nonsingular. In fact, with the definitions
A derivation of this result is indicated in Problem 2.4.2. Now the coefficient matrix in Equation (2.4.11) can be written as
and hence it is nonsingular if
Using (2.4.12) and (2.4.7) it follows that
where
Furthermore, just as in the proof of Theorem 1.2.6, we find from (2.4.10b) and the form of W that
where k(x) is, as usual, the bound
We see that taking all the ∆j “small” reduces the bound in (2.4.14b). But since 
, this implies “large” values for J which increases the bound in (2.4.14a). If we assume that k(x) is bounded on [a, b] and take subintervals of equal length, ∆j= (b – a)/J, j= 1, 2, …, J, then the in-equality (2.4.13) is satisfied if
On the other hand, Theorem 1.2.6 is valid if we replace the hypothesis (1.2.18e) by the stronger condition
For large J and m the above inequalities are almost equivalent, since then
We should also observe that when (2.4.13) is satisfied, the proof of Theorem 1.2.6 applies to the boundary-value problem (2.4.5) to show that it has a unique solution.
What we have shown above is, essentially, that the parallel shooting procedure does not make matters worse than they were for the original initial-value problem. We have not been able to show, for instance, that the length of the interval [a, b] over which our proof of existence and uniqueness holds, can be increased. But the purpose of the parallel shooting method is to overcome certain practical difficulties that occur in ordinary shooting. While actual computations verify its usefulness in special cases, we may also give some theoretical justification as follows. In solving the initial-value problems (2.3.2) and (2.4.8) numerically, errors are of course introduced. If the same stable scheme of order p, say, is used in each case with the same net spacing h, and the functions f(x, u) and F(t, U) are sufficiently smooth, then the numerical errors at points of the intervals [a, b] and [0, 1] can be bounded by
Some estimates of this form are derived in Section 1.3. The constants M1 and M2 depend in a simple algebraic manner on bounds for the functions f and F and possibly some of their derivatives. However, the constants K1 and K2 are essentially the Lipschitz constants for f and F, respectively, or polynomials in these constants. Then if we take ∆j = |b – a | J, j = 1, 2, …, J, in the subdivision (2.4.2) it is clear from (2.4.3) and (2.4.6) that
Then since solutions of (2.3.2) at x = b correspond to those of (2.4.8) at t = 1, we see that the bound on the error in the numerical solution is reduced exponentially by parallel shooting [that is, the bound is proportional to (exp [K1|b – a|])1/J rather than to (exp[K1|b – a|])]. The constant M2 is at worst larger than M1 by a factor which is a small power of J. To show the detailed dependence of the quantities M1, M2, K1, and K2 on bounds for f, F, and their derivatives is rather complicated. Problems 1.3.3–6 contain some results which should clarify matters. More details can be found in Isaacson and Keller, Chapter 8, Sections 2 and 3. Very explicit error estimates are contained in Henrici (1962) and (1963).
A schematic diagram of the parallel shooting scheme presented above is shown in Figure 2.4.1. By analogy it is easy to indicate other parallel shooting schemes which may in many cases be more efficient. For instance, if we agree to take J = 2N, so that there is always an even number of subintervals, we can shoot from each odd-labeled subdivision point toward each of its two adjacent even-labeled points. A glance at Figure 2.4.2 for the special case J =4should suffice. Here there are only two n-vectors s1 and s2 to be determined [by applying continuity at x2, that is u2(1) = u3(1), and the appropriate boundary conditions on u1(1) and u4(1)]. The continuity conditions at x1 and x3are automatically satisfied. We do not bother with a detailed formulation of this case as the diagram should suffice. In this manner the systems replacing (2.4.9) and (2.4.11) are only half the order they would be if all shooting were in the same direction for the same number of subintervals.
Figure 2.4.1
Figure 2.4.2
A rather common special case of the above modification occurs for J=2 and is generally described in terms of shooting from each endpoint and matching at some interior point. In fact all of our schemes are just as well considered to be multiple matching generalizations of this well-known case.
If it is known that all solutions of a given system of differential equations decay in one direction then all integrations should be done in that direction. However if solutions grow more or less equally in both directions then there is no preferred direction and the more efficient process of shooting in both directions should be used.
A rough practical guide to the choice of subintervals to be used in parallel shooting can be based on the old-fashioned technique (when calculating with fixed point arithmetic) of rescaling the solution as it grows. That is, integrate the initial-value problem
using some “typical” initial value ξ on some sufficiently fine net {xv} with x0 = a. After each step of the integration test to see if the growth condition
where R is some selected reasonable growth factor, say R = 103 (the choice of R depends upon the machine word length, as well as other factors). When this test is passed, take xv to be one of the points in (2.4.2a), replace ||y(x0)|| in the test by the new value, ||y(xv)||, and continue the integration and testing. Some experimenting with choices of ξ and R may be required in very difficult cases.
Power series are occasionally applied or suggested as a means for solving boundary-value problems. Their application is in fact a variant of the initial-value method. For example, if f(x, y) in Equation (2.4.1a) is an analytic function of the n + 1 variables (x, y), then we may seek power-series solutions in the form
Using this in the expanded form of Equation (2.4.1a) and equating coefficients of like powers of (x–a), we obtain (recursive) formulae for the coefficients ak, k = 1, 2, …, in terms of the initial point a, and the leading coefficients a0 [i.e., ak = ak(a0)].
Now a0 is to be determined so that the boundary conditions (2.4.1b) are satisfied. This yields the system
Since the ak for k ≥ 1 are functions of a0 the system (2.4.16) can be considered a transcendental system of equations in a0. [Of course if there are p independent conditions on y(a) in (2.4.1b), or equivalently on a0 in (2.4.1b), then the system is easily reduced to one containing only q = n – p equations and unknowns.]
However, if power series are employed to solve the initial-value problem (2.3.2) then the resulting transcendental system ϕ(s) = 0, with ϕ(s) defined in (2.3.3), is identical to the system (2.4.16) if we set s ≡ a0. In this sense the power-series procedure is but a special case of the initial-value method.
To justify the power-series method in any case requires a study of the convergence of the appropriate series (2.4.15) over the interval a ≤ x ≤ b. Thus, at best, the method has rather limited applicability. In actual calculations we do not use (2.4.15) but some truncated version, say,
The error committed in this way is formally analogous to the truncation error introduced when solving the initial-value problem (2.4.1) by some numerical scheme. While the calculations using (2.4.17) are always well-defined, the formal result as N → ∞ may be nonsense if the radius of convergence is too small, that is, less than b – a. (As indicated in Section 2.3 the schemes based on direct numerical integration of the initial-value problem have no such limitation imposed by radius-of-convergence difficulties.)
However, there are several ways in which the power-series method may be salvaged. One way is by direct analytic continuation from a to b by expansions about a sequence of intermediate points, xj, j = 1, 2, …, J. This procedure is essentially the same as using the one-step integration scheme known as the Taylor expansion method to solve the initial-value problem (2.3.2). When the motivation is power series one would employ many terms in the expansions of form (2.4.17) and few “expansion points” xj, that is, small J. Conversely, many closely-spaced “net points,” xj, and few terms in the Taylor series would be used in schemes motivated by numerical integration. The difference in point of view depends, of course, on which limiting process is intended, maxj |xj, – xj+1| → 0, or N → ∞; only in the latter case is analyticity required.
Some other ways in which power series may be used are suggested by the parallel shooting variants discussed above. For instance, the solution can be expanded about each endpoint and the two series matched at some common interior point. Then two unknown leading-coefficient vectors, say a0 and b0, must be determined by the matching conditions and the boundary conditions. Obviously the solution could be expanded about N interior points, matched at N – 1 interior points separating these expansion points, and the boundary conditions imposed to determine N unknown leading-coefficient vectors. For instance, the odd-labeled points in Figure 2.4.2 may be taken as the origins for expansions, and at the even-labeled interior points the appropriate series would be matched. The convergence problems and questions of accuracy in such power-series methods are untreated in the literature.
2.4.1 The problem y″ = m2 sinh q2y, y(0) = y(1) = 0 has the unique solution y(x) ≡ 0. Show that the initial-value problem
has a singularity located at approximately
Thus in order that x∞ > 1, that is, lie outside [0, 1], the initial parameter s must satisfy
For m = q = 4 the initial parameter s must be accurate to within e–16. [This problem is from B. A. Troesch (1960)], [HINT: By integrating 
, the solution can be obtained in the implicit form: 
. Setting u = ∞, we obtain an expression for x∞ which can be estimated by expanding the integrand appropriately.]
2.4.2 For the matrix (P + Q) determined by Equation (2.4.7) verify the block factorization:
Use the above to deduce the result in (2.4.12).
2.4.3 Using the fact that W = diag {Wj} and the assumption that the nth-order matrix (A + BWJ ··· W1) is nonsingular, derive an explicit representation for the inverse of the nJ-order matrix (P + QW). [HINT: Factor (P + QW) into block-triangular form with first row (A + BWj ··· W1), BWj ··· W2,···,BWJ in a manner similar to that of Problem 2.4.2. Find the inverse as the product of triangular matrices in analogy with (2.4.12).]
2.4.4 Formulate the parallel shooting scheme implied by Figure 2.4.2 (generalized to J = 2N) for the boundary-value problem (2.4.1). Examine the linear system obtained by applying Newton’s method to solve for the nN-dimensional vector s.
2.4.5 Consider the boundary-value problem of order n = 2m
where y, z, f, g, α, and β are m-dimensional and the mth-order matrices A and B are nonsingular. Formulate the parallel shooting schemes of both types, Figures 2.4.1 and 2.4.2, for the same J = 2N intervals for this problem. Compare the application of Newton’s method in each case and formally determine the relevant inverse matrices. Discuss the solution when f andgare linear in yand z.
Section 2.1 A technique for reducing the effect of cancellation errors in some linear boundary-value problems is given in Godunov (1961). This procedure has been made more practical by S. Conti (1966), whose scheme can be formulated as a special form of parallel shooting.
Section 2.2 Shooting methods for nonlinear second-order problems are described by Collatz (1960) and Fox (1957).
Section 2.3 The Newton method applied to shooting with somewhat specialized boundary conditions is described in Goodman and Lance (1956). Rather than solve a small-order linear system to obtain the corrections to the initial data they integrate (numerically) the adjoint system backwards.
Since the variational systems are linear, predictor-corrector schemes are better replaced by the corrector alone, properly solved [see Riley, Bennett, and McCormick (1967)].
Section 2.4 A limiting form of parallel shooting yields a finite-difference scheme (that is, when the modified Euler one-step method is used over each interval). Multiple shooting in various forms has been considered by, among others: Fox (1960), Morrison, Riley, and Zancanaro (1962), and Kalnins (1964). For linear problems, Conti’s (1966) procedure is closely related. However, one of the most important features of parallel shooting— the exponential reduction in error growth—does not seem to have been recognized previously.
For an example of the use of power series in solving nonlinear boundary-value problems see Reiss, Greenberg, and Keller (1957) or Weinitishke (1960), who employs analytic continuation.
