  We present a straightforward and reliable continuous method for computing the
full or a partial Lyapunov spectrum associated with a dynamical system
specified by a set of differential equations. We do this by introducing a
stability parameter beta>0 and augmenting the dynamical system with an
orthonormal k-dimensional frame and a Lyapunov vector such that the frame is
continuously Gram-Schmidt orthonormalized and at most linear growth of the
dynamical variables is involved. We prove that the method is strongly stable
when beta > -lambda_k where lambda_k is the k'th Lyapunov exponent in
descending order and we show through examples how the method is implemented. It
extends many previous results.
