I have some a function that I need to solve numerically over an entire domain (Ideally I would like an InterpolatingFunction object as output), which does not include any derivatives. I have found a way to do it using NDSolve, however I feel like there should be a better solution. Here is a simple example to show my NDSolve method:

Solving:

$y^2+y=x$

The solutions:

$y(x) = frac{1}{2}left(-1 pm sqrt{1+4x}right)$

Numerically Solving using NDSolve:

`NDSolve(D(y(x)^2 + y(x), x) == D(x, x) && y(0)^2 + y(0) == 0, y(x), {x, -5, 5}) //Quiet`

This operation gives me two InterpolatingFunctions that, when plotted, match the analytical solutions, but feels very hacky. I have also noticed that for some functions it does not produce solutions for the entire function domain. Is there a better way to do this?