This is an amplification of my comment on Vladimir's answer. It's actually not at all hard to see the generality of the surface metrics that admit a $k$-th degree polynomial first integral of their geodesic flow. Here is a sketch:
The question is local, so we can look at metrics on an open set $U\subset\mathbb{R}^2$. Moreover, modulo diffeomorphism, we can assume that the metric is conformal, i.e., $g = e^{2u}\bigl(dx^2+dy^2\bigr)=e^{2u}dz\circ d\bar{z}$, with cometric
$$
\hat g = e^{-2u}\,\frac{\partial}{\partial z}\circ\frac{\partial}{\partial \bar{z}},
$$
which is a function on the symplectic manifold $T^*\mathbb{R}$.
Now, a polynomial first integral of degree $k$ is a function $p$ on the tangent bundle of $\mathbb{R}^2$ of the form
$$
p = v_0(x,y)\,dx^k + v_1(x,y)\,dx^{k-1}dy + \cdots + v_k(x,y)\,dy^k.
$$
Let $\hat p:T^*\mathbb{R}^2\to\mathbb{R}$ be its $g$-dual, considered as a function on $T^*\mathbb{R}$.
The condition that $p$ be constant on the geodesic flow of $g$ is simply that $\hat g$ and $\hat p$ Poisson commute, i.e., $$ \left\{\hat g, \hat p\right\} = 0. $$
Since the expression $\left\{\hat g, \hat p\right\}$ is polynomial of degree $k{+}1$ in the momenta, this is $k{+}2$ first-order equations for the $k{+}2$ unknowns $u, v_0,\ldots, v_k$. It is not difficult to see that this system can be placed in Cauchy-Kowalewskaya form, so analytic solutions are determined by specifying these $k{+}2$ functions analytically along an appropriately non-characteristic curve.
However, there is still too much symmetry in the problem, namely the conformal transformations of the complex plane. Generically, you can get rid of this as follows. If you write $\hat p$ in the form $$ \hat p = h_0(z,\bar z)\ \left(\frac{\partial}{\partial z}\right)^k + h_1(z,\bar z)\ \left(\frac{\partial}{\partial z}\right)^{k-1} \circ \frac{\partial}{\partial {\bar{z}}} + \cdots + h_k(z,\bar z)\ \left(\frac{\partial}{\partial {\bar{z}}}\right)^k $$ where $\overline{h_j} = h_{k-j}$, you find that the vanishing of the Poisson bracket implies that $h_0$ is actually holomorphic. You can always assume that $h_0$ is not identically vanishing because, otherwise, you could factor out a copy of $\hat g$ from $\hat p$ and so reduce the order of the integral. Now, use a holomorphic transformation to make $h_0\equiv1$. This reduces the number of unknowns by $2$ and reduces the number of equations by $2$, and the resulting first order system of $k$ equations for $k$ unknowns now has exactly the right generality, showing that the general solution depends on $k$ functions of one variable.
All of what I have written was known more than a century ago, and you can find an account of it in Volume 3 of Darboux', Leçons sur la Théorie Générale des Surfaces et les Applications Géométriques du Calcul Infinitésimal. Possibly, what Vladimir is referring to about the difficulties is something else: Because we now know that there are restrictions on a sufficiently high order jet of a metric in order for it to admit a geodesic first integral of order $k$, we might want to know what those conditions are, explicitly in terms of some known invariants. However, this turns out to be extremely complicated once you get beyond degrees $1$ and $2$.