We want to solve 
. Discritize the
interval [a,b] into N+1-subintervals defined by 
. At each mesh point, approximate derivatives by finite differences,
substitute into the ODE, and solve the resulting nonlinear system of equations
for the unknowns 
.
It appears this method is easy. Of course we find ourselves confronted with all
the difficulties we encountered earlier with nonlinear systems. It is often helpful
to parameterize the ODE with some parameter, say 
, so that 
corresponds to an equation for which you know the exact solution (or its real easy to
compute), and 
corresponds to the desired ODE. Then one can slowly increase
the value of the parameter, solving the parameterized system along the way.
Note that if you use an iterative solver, you don't have to be too accurate
in your early approximate solutions; you can crank down the tolerance as the
parameter iteration proceeds.