In two recent papers [1, 2], we have developed an algorithm for the solution of linear differential systems of the type

formula here

on an interval [X, ∞). Here Z is an n-component vector, ρ is a scalar factor, D is a constant diagonal matrix and R is a perturbation matrix such that R(x)=O(x−δ) as x→∞, for some δ>0.The algorithm implements a repeated transformation process by means of which (1·1) is transformed into other systems whose perturbation matrices are of successively smaller orders of magnitude as x→∞. When the perturbation reaches a prescribed accuracy, the Levinson asymptotic theorem [4, 5] provides the solution of the final system in the process, the solution also possessing that accuracy. The corresponding solution of (1·1) is then obtained by transforming back.

The algorithm has two main aspects. The first is the algebraic one of generating symbolically the matrix terms which appear in the transformation process. In [1] and [2], the algorithm was set up in such a way that this aspect is implemented in the symbolic algebra system Mathematica. The other aspect is the numerical one of deriving from the algebra the values of the solutions of (1·1), given ρ, D and R. Included here are the determination of explicit error-bounds and procedures for handling large numbers of matrix terms.

In this paper we develop these ideas in two different but not unrelated directions. First, we replace the constant matrix D by one which has entries of differing orders of magnitude as x→∞. Our methods are sufficiently indicated by two orders of magnitude, and thus we consider the system

formula here

where formula here

with constant and diagonal D˜ and D, while λ is a scalar factor such that λ(x)→∞ as x→∞. The amplification of our previous algorithm to cover this new situation is dealt with in Sections 2–4 below. At the end of Section 3, we state the Levinson asymptotic theorem in the form that we need, with the necessary accuracy incorporated. As in [1] and [2], a typical ρ(x) is xa (α>−1). However, the value α=−1 also occurs in applications, and our second extension is to include α=−1 in the scope of the algorithm. This matter is the subject of Section 5.