This paper is concerned with boundary value problems for singularly perturbed linear systems of the form \(\epsilon y'=A_{1,1}y+A_{1,2}z+f_ 1,\) \(z'=A_{2,1}y+A_{2,2}z+f_ 2,\) \(B_ 0(\epsilon)\left( \begin{matrix} y\\ z\end{matrix} \right)(0)+B_ 1(\epsilon)\left( \begin{matrix} y\\ z\end{matrix} \right)(1)=\beta.\) \(A_{1,1}(t,\epsilon)\) is assumed to satisfy \(E(t)^{- 1}A_{1,1}(t,0)E(t)=diag(\lambda_ 1(t),...,\lambda_ n(t))\) where the \(\lambda_ i(t)\) stay away from the imaginary axis all over [0,1]. For the numerical solution of this problem, collocation methods are discussed which are based on Gauss or Lobatto type Runge-Kutta formulae. These methods have been treated in Part I of this work [SIAM J. Numer. Anal. 20, 537-577 (1983; Zbl 0523.65064)] for the case of singular problems with constant coefficients, whereas the second author [Math. Comp. 28, 449-464 (1974; Zbl 0284.65067)] is especially concerned with the trapezoidal and box schemes for the numerical solution of the more general problem. - The present paper gives some information on the implementation of both the Gauss and Lobatto schemes. Since the differential equation has layers at both boundaries, for the numerical integration a finer mesh is needed near the boundaries than in the interior of the interval. The authors give explicit formulas for the meshes. The major part of this paper is taken up by the error analysis. For this, both the solutions to the differential equation and to the discrete problem are decomposed into components each showing a typical asymptotic behaviour. Moreover, due to the layers, this is done on three subintervals separately.