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$.

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 quasilinear first order system can locally be placed in Cauchy-Kowalewskaya form, so analytic solutions are determined by specifying these $k{+}2$ functions analytically along an appropriately non-characteristic curve. (Not all solutions are real-analytic, however.)

Now, there is still too much symmetry in this formulation, namely the conformal transformations of the complex plane. Generically, one can get rid of this as follows. If one writes $\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}$, one finds that the vanishing of the Poisson bracket implies that $h_0$ is actually holomorphic.

Now, one can always assume that $h_0$ is not identically vanishing because, otherwise, one could factor out a number of copies 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$ (as it fixes the real and imaginary parts of $h_0$) and reduces the number of equations by $2$ (because two of the equations are now identities), and the resulting first order system of $k$ equations for $k$ unknowns now has exactly the right generality and can still be put in C-K form locally, showing that the general (analytic) solution depends on $k$ functions of one variable.

All of what I have written was known more than a century ago, and one 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 polynomial geodesic first integral of degree $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 one goes beyond degrees $1$ and $2$.

What one can say is this (which gives a more precise version of the Theorem that Vladimir states in his answer): Consider the space $\mathsf{G}_\ell$ consisting of $\ell$-jets of surface metrics modulo diffeomorphism. This is a (singular) space that is finite dimensional, of dimension $1$ when $\ell=2$ and of dimension $\tfrac12(\ell{-}2)(\ell{+}1)$ when $\ell\ge 3$. (The singular locus of $\mathsf{G}_\ell$ has properly smaller dimension, and is nonempty for $\ell\ge 3$.) Carefully applying the above analysis shows that, when $\ell\ge2k{+}2$, the locus $\mathsf{F}_\ell(k)\subset \mathsf{G}_\ell$ of diffeomorphism classes of $\ell$-jets of metrics that admit a nontrivial polynomial geodesic first integral of degree $k$ or less has codimension $C_\ell(k) = \tfrac12(\ell{-}2k{-}2)(\ell{+}1)+1$ in $\mathsf{G}_\ell$.

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. $$

Now, the expression $\left\{\hat g, \hat p\right\}$ is always polynomial of degree $k{+}1$ in the momenta, so 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.

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$.