  We consider the least-squares finite element method (lsfem) for systems of
nonlinear ordinary differential equations and establish an optimal error
estimate for this method when piecewise linear elements are used. The main
assumptions are that the vector field is sufficiently smooth and that the local
Lipschitz constant, as well as the operator norm of the Jacobian matrix
associated with the nonlinearity, are sufficiently small when restricted to a
suitable neighborhood of the true solution for the considered initial value
problem. This theoretic optimality is further illustrated numerically, along
with evidence of possible extension to higher-order basis elements. Examples
are also presented to show the advantages of lsfem compared with finite
difference methods in various scenarios. Suitable modifications for adaptive
time-stepping are discussed as well.
