  We propose an algebraic method that finds a sequence of functions that
exponentially approach the solution of any second-order ordinary differential
equation (ODE) with any boundary conditions. We define an extended ODE (eODE)
composed of a linear generic differential operator that depends on free
parameters, $p$, plus an $\epsilon$ perturbation formed by the original ODE
minus the same linear term. After the eODE's formal $\epsilon$ expansion of the
solution, we can solve order by order a hierarchy of linear ODEs, and we get a
sequence of functions $y_n(x;\epsilon,p)$ where $n$ indicates the number of
terms that we keep in the $\epsilon$-expansion. We fix the parameters to the
optimal values $p^*(n)$ by minimizing a distance function of $y_n$ to the ODE's
solution, $y$, over a given $x$-interval. We see that the eODE's perturbative
solution converges exponentially fast in $n$ to the ODE solution when
$\epsilon=1$: $\vert y_n(x;\epsilon=1,p^*(n))-y(x)\vert<C\delta^{n+1}$ with
$\delta<1$. The method permits knowing the number of solutions for Boundary
Value Problems just by looking at the number of minima of the distance function
at each order in $n$, $p^{*,\alpha}(n)$, where each $\alpha$ defines a sequence
of functions $y_n$ that converges to one of the ODE's solutions. We present the
method by its application to several cases where we discuss its properties,
benefits and shortcomings, and some practical algorithmic improvements.
