Systematic way of obtaining conservation laws in dynamical systems Motivation
Consider a point particle of mass $m$ moving in $\mathbb{R}^3$ under the influence of some force field $\vec{F}(\vec{r},t)$. The fundamental equation governing the dynamics of this system is Newton's second law stating $$\boxed{\vec{F} = \frac{d \vec{p}}{dt}}$$ where $\vec{p} = m \vec{v}$ is the (linear) momentum of the particle. This is also an example of a conservation law. It states that a free particle, i.e. one for which there is no force acting on it, moves in a way such that its (linear) momentum is conserved $$\vec{F} = \vec{0} \quad \iff \quad \vec{p} = constant.$$ In this question, I am interested in a systematic way of constructing all such conserved quantities.
My approach
My approach, in a nutshell, is to multiply the Newton's law with various quantities and, after some manipulation, recognize that some quantities are constant when $\vec{F} = \vec{0}.$ In the following, I will list some examples.


*

*Dot product $\vec{v} \cdot \vec{F}$ leads to the conservation of kinetic energy $T = \frac{1}{2}m \vec{v}^2$. We have $$\boxed{\vec{v} \cdot \vec{F} = \frac{d}{dt} \left( \frac{1}{2}m \vec{v}^2 \right)}$$ so that $$\vec{F} = \vec{0} \quad \iff \quad T = constant.$$

*Cross product $\vec{v} \times \vec{F}$ apparently does not lead to any conservation laws because we trivially have $$\boxed{\vec{v} \times \vec{F} = m \vec{v} \times \vec{a}}$$ and the RHS cannot be written as a time derivative. Therefore, it looks as though $\vec{v} \times \vec{F}$ is not a useful quantity. (I could be wrong about this, so please correct me.)

*Dot product $\vec{r} \cdot \vec{F}$ does not lead to a conservation law due to an extra term on the RHS, $$\boxed{\vec{r} \cdot \vec{F} = \frac{d}{d t}(\vec{r} \cdot \vec{p})+2T}$$ However, this equation is the starting point in the derivation of Virial theorem, so it's not totally useless.

*Cross product $\vec{r} \times \vec{F}$ leads to the conservation of angular momentum $\vec{L} = \vec{r} \times \vec{p}$ because $$\boxed{\vec{r} \times \vec{F} = \frac{d}{d t} (m \vec{r} \times \vec{v})}$$ so that $$\vec{F} = \vec{0} \quad \iff \quad \vec{L} = constant.$$

*Multiplication with time leads to the conservation of the center of mass motion, $$\boxed{t \vec{F} = \frac{d}{dt} (t \vec{p} - m \vec{r})}$$ so that $$\vec{F} = \vec{0} \quad \iff \quad t \vec{p} - m \vec{r} = constant.$$
The question
In my approach, I have been able to reproduce the 10 well known integrals of motion simply by guesswork. However, I have no guarantee that there aren't any more integrals of motion. So, is there a systematic way of finding all integrals of motion for a given dynamical system? 
 A: First off, read up on Noether's theorem and especially the Hamiltonian and Poisson formulations of it to get up to speed with modern tools. In a nutshell, symmetries and integrals-of-motion exist in a 1-to-1 correspondence, although things get hairy when there's dependence involved in which case you choose isotropy subgroups of symmetries. I've personally found identifying symmetries and then using Noether's theorem to find the integrals-of-motion to be easiest strategy instead of directly searching for the integrals themselves.
As far as I know, there is no known method for determining the total number of integrals-of-motion or symmetries a given system may possess. This is because integrals-of-motion and symmetries are characteristics of global system behavior, not just local behavior. To find or somehow guarantee the existence of symmetries and integrals-of-motion, you would have to also know all of the global behavior, which there are plenty of examples where this is impossible.
Your examples are more or less a guess and check approach, which is in fact one of the best methods available for this kind of work. This is a numerical approach that you might find as a satisfactory answer. Calculus of variations and optimal control are just generalizations of the mechanics you are used to. I've used this technique before and have had mixed results, such as it will find 1 of 3 symmetries and integrals-of-motion that I already know exist. The fact that such a numerical procedure exists tells you that this is a tough problem with no known general answer, and the fact that the numerical procedure only works so-so tells you this is an area we need some more research done.
Edit: I should also mention that symmetries are not the same as infinitesimal symmetries. If I had to speculate that such a process existed for determining all symmetries a system may possess, it would be done at the level of the infinitesimal symmetries, which are a Lie algebra of vector fields. Then Lie III would identify the full group of symmetries. This is pure speculation on my part and I no evidence to back it up.
