Did you try softwares doing rigorous numerical computation? For example arb is a C library freely available in Sage. Such software gives guaranteed enclosures $f(I) \subset J$ where $I,J$ are intervals. So an obvious (but not necessarily optimal) algorithm for finding the maximum of your function $f:I \to \mathbb{R}$ is to subdivide your interval $I$ and apply $f$ to each subinterval, and keep subdividing until you get the maximal value to the precision you want. Personally I would judge this as rigorous and acceptable in a proof. This should work as long as the functions you need are present in arb, at worst the algorithm will be too slow.