The standard shooting and fitting algorithm for non-linear two-point boundary value problems derives from conventional coordinate perturbation theory. We generalize the algorithm using the renormalized perturbation theory of strained coordinates. This allows for the introduction of an arbitrary function, which may be chosen to improve numerical convergence. An application to a problem in stellar structure exemplifies the algorithm and shows that, when used in conjunction with the standard procedure, it has superior convergence compared to the standard one alone.