The framework of sum ofsquares / semidefinite programming (SOS / SDP) allows one to compute arbitrarily good approximations to such problems. The general idea is as follows.
First, we write the problem as an optimization problem whose decision variables are the polynomial's coefficients and some other auxiliary variables. Second, we transform this into a conic program with linear objective function, linear equality constraints, and some conic constraints that certain decision variables define the coefficients of a polynomial which is everywhere nonnegative.
If we ignore for a moment the issue of nonnegativity on an interval / box, your problem would become something like
\begin{align}
minimize & \epsilon \\
\text{subject to} & P-F \text{ is a nonnegative polynomial and } \\
& \epsilon - P \text{ is a nonnegative polynomial.}
\end{align} Note that, for example, $\epsilon$ is a decision variable so by introducing new variables and equality constraints, we can implement the second line with the types of constraints mentioned above by saying that $Q = \epsilon - P$ (linear equality constraint on coefficients) and $Q$ is nonnegative (nonnegative cone constraint).
Third, we use the fact that while the cone of nonnegative polynomials is not analytically tractable (except for certain small numbers of variables or degrees), it has nice inner approximations. The simplest are the sums of squares: a polynomial which is a sum of squares of other polynomials (hereafter a "sum of squares") is always nonnegative. Replacing the cone of nonnegative polynomials everywhere by the cone of sums of squares, we get an inner approximation to the desired optimization problem. That is to say, the optimal solution will be feasible for the original problem and the corresponding $\epsilon$ will be an upper bound on the optimal $\epsilon$ for the original problem.
Fourth, there is a purely syntactic procedure to convert such a "sum of squares program" to a semidefinite program (i.e. an optimization problem with symmetric matrix decision variables, linear objective and equality constraints, and positive semi definiteness constraints). This can then be solved by an SDP solver such as SeDuMi or SDPT3, and the result translated back to give the coefficients of the optimal polynomial.
That is the general scheme. Similar ideas work for constraints that polynomials be nonnegative on intervals, boxes, etc. For example, a univariate polynomial $p$ is nonnegative on $[-1,1]$ if and only if there exist sums of squares $s,t$ such that $p(x) \equiv s(x) + (1-x^2)t(x)$.
Also, there are tighter relaxations than the sums of squares if the approximations end up being too loose. For a lot more detail on all of this, see these lecture notes, especially Lecture 10.
Fortunately there is software to take care of all of this automatically in MATLAB (perhaps I should have mentioned this first…). I recommend googling SOSTOOLS and YALMIP for a start. These are well-developed enough that you do not need to know all the theory in order to use them effectively, but it is worth at least knowing what I wrote above to know that the answers will sometimes be approximations.
