  Four generalizations of the Phase Integral Approximation (PIA) to sets of N
ordinary differential equations of the Schroedinger type: u_j''(x) + Sum{k = 1
to N} R_{jk}(x) u_k(x) = 0, j = 1 to N, are described. The recurrence relations
for higher order corrections are given in the form valid in arbitrary order and
for the matrix R_{jk} either hermitian or non-hermitian. For hermitian and
negative definite R matrices, the Wronskian conserving PIA theory is formulated
which generalizes Fulling's current conserving theory pertinent to positive
definite R matrices. The idea of a modification of the PIA, well known for one
equation: u''(x) + R(x) u(x) = 0, is generalized to sets. A simplification of
Wronskian or current conserving theories is proposed which in each order
eliminates one integration from the formulas for higher order corrections. If
the PIA is generated by a non-degenerate eigenvalue of the R matrix, the
eliminated integration is the only one present. In that case, the simplified
theory becomes fully algorithmic and is generalized to non-hermitian R
matrices. General theory is illustrated by a few examples generated
automatically by using author's program in Mathematica, published in
arXiv:0710.5406.
