The simple answer, overall, is that you always have to assume an integral of motion in a certain form and look if the conditions for its existence are even fulfilled. Generically, they are not, and if no explicit symmetry indicates full integrability, one cannot expect that the motion will be fully integrable; on the contrary, parts of the phase space will almost always be chaotic. Methods of finding the "hidden" integral of motion vary on a case to case basis and I will discuss only the case of looking for the "last" integral of motion needed for full integrability given that we already have integrability coming from stationarity and axisymmetry.
Additionally, I will interpret your question in two ways, 1) the motion of a relativistic particle, and 2) the motion of a Newtonian particle in a potential field. I will also discuss only global integrals of motion, i.e. those valid for any trajectory.
As for 1), you have to satisfy four-velocity normalization $g^{\mu \nu} u_\mu u_\nu = -1$. However, the Hamiltonian you give leads to a normalization $g^{\mu \nu} p_\mu p_\nu = -\mu^2 - V(q)$ with $\mu$ some constant which means that (on-shell), your momenta are related to four-velocity as $p_\mu = \sqrt{\mu^2 + V} \,u_\mu$.
Further computations show you that your Hamiltonian in fact reproduces trajectories parametrized by a non-affine parameter $\lambda$ related to proper time $\tau$ as $d \tau =\sqrt{\mu^2 + V} d\lambda $. The Hamiltonian which reproduces trajectories parametrized by proper time would be
$$H = \frac{1}{2 \sqrt{\mu^2 + V}} g^{\mu \nu} p_\mu p_\nu - \frac{1}{2} \sqrt{\mu^2 + V} $$
However, we can push further, and reparametrize by another parameter $d \tau = dp/\sqrt{\mu^2 + V} $ which leads us to the Hamiltonian
$$H = \frac{1}{2 (\mu^2 + V)} g^{\mu \nu} p_\mu p_\nu $$
I.e., what you are looking for is the integral of motion for geodesics in the metric $\tilde{g}^{\mu \nu} = g^{\mu \nu}/(\mu^2 + V)$. For that, you must look for Killing tensors of the metric $\tilde{g}^{\mu \nu}$, as discussed in this answer.
Either way, it should be understood that the only four-dimensional vacuum solution to Einstein equations possessing an irreducible (not trivially composed from simple Killing vectors) Killing tensor is the Kerr-NUT-(A)dS space-time. If you have a generic stationary axisymmetric vacuum metric, the motion will simply not be integrable.
As for 2), you are also allowed to do this reparametrization trick, but your metric will be valid only on isoenergetic surfaces in phase space because do not have anything like four-velocity normalization in Newtonian physics. The metric obtained by this method is called the Jacobi metric (you can find more in Pettini's Geometry and Topology in Hamiltonian Dynamics) and you may look for the Killing tensor of this metric to find the other integral of motion.
However, it is more convenient to assume an integral of motion of the form
$C = K^{ij}p_i p_j + K(q)$
and look for consequences of the requirement that $\dot{C} = 0$ in Cartesian coordinates without any funny reparametrizations etc. For your Hamiltonian and in Cartesian coordinates you get (working only with $i,j$ indices to indicate the Newtonian nature of the problem) the conditions
$$\partial^{(k}K^{ij)} = 0$$
$$\partial^{j} K = 2 K^{jk} \partial_k V$$
Where we now have just partial instead of covariant derivatives which makes it possible to obtain explicit conditions on $V$ to generate integrable motion. When the dust settles, you obtain that the motion in a stationary and axisymmetric $V$ will be integrable if and only if $V$ is of the form
$$V = \frac{C_1}{\sqrt{x^2 + y^2 + (z-C_2)^2}} + \frac{C_1}{\sqrt{x^2 + y^2 + (z+C_2)^2}} + C_3 (x^2 + y^2 + z^2)$$
where $C_1,C_2,C_3$ are some constants. If you have an axisymmetric stationary potential, it can be shown (with quite a lot of work especially for the quartic case) that a polynomial integral of motion of order higher than quadratic cannot exist.
I.e., you have to take a look at your axisymmetric stationary potential and see whether it looks like the one given above. If not, it will not be globally integrable.
This all is nicely summarized in the paper by Charalampos Markakis from 2014, you will also find good references in that paper.
