Four generalizations of the phase integral approximation (PIA) to sets of ordinary differential equations of Schrödinger type [uj″(x)+∑k=1NRjk(x)uk(x)=0, j=1,2,…,N] are described. The recurrence relations for higher order corrections are given in a form valid to arbitrary order and for the matrix R(x)[≡{Rjk(x)}] either Hermitian or non-Hermitian. For Hermitian and negative definite R(x) matrices, a Wronskian conserving PIA theory is formulated, which generalizes Fulling’s current conserving theory pertinent to positive definite R(x) matrices. The idea of a modification of the PIA, which is 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 nondegenerate eigenvalue of the R(x) 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(x) matrices. The general theory is illustrated by a few examples automatically generated by using the author’s program in MATHEMATICA published in e-print arXiv:0710.5406 [math-ph].