Let $f:\mathbb{R}_0^+\to \mathbb{R}$ be defined by some combination of the four basic operations and square roots. (The argument of square-roots is assumed is to be non-negative, and the value of square roots is defined to be non-negative as well.) How can one find the global maximum and minimum of $f$ numerically, with full rigor?

(A rigorous numerical solution has to be given as an interval, e.g., [1.500001,1.500002], within which the maximum (say) lies.)

Rather convincing but not quite rigorous method:

Plot $f'(x)$, see roughly where the zeroes are, and narrow them down using some standard (e.g. SAGE, mathematica). Then verify this (non-rigorous) datum rigorously simply by checking that $f'(x)$ has different signs at $x=x_0-\epsilon$ and at $x=x_0+\epsilon$, where $x_0$ is an alleged zero of $f'(x)$. Compare the values of $f(x)$ at all $x_0$, and also at $x=0$ and as $x\to \infty$.

There is really only one way in which this procedure isn't rigorous: there is no guarantee that there aren't any zeroes of $f'(x)$ we have missed. (If $f$ is a polynomial, this can be easily dealt with, but $f$ is not necessarily a polynomial.)

Are there any (free) programs out there that for the minimum and maximum of $f$ rigorously?