In Hamiltonian mechanics we search change of coordinates that leaves the Hamilton equation invariant: these are the canonical transformations.
My question is: why we want to leave the equations invariant? I mean: we want to solve a differential equation (the Hamilton equation), so why restrict ourselves only to a certain kind of change of coordinates? Why we care about the "form"? I suppose that sometimes could exist some non-canonical transformation that makes me solve the equation... or not?
See What's the point of Hamiltonian mechanics? for a closely related question. The answers in that thread mention various benefits with the Hamiltonian formulation, among other things:
Analysis of structural patterns, symmetries & conservation laws, separability & integrability, Liouville theorem, Hamilton-Jacobi theory, etc, without solving the EOMs.
Contact to (mainly non-relativistic) quantum mechanics and statistical physics.
If we are only interested in solving the EOMs, then the Lagrangian formulation is usually easier.
In this answer, we assume that by a canonical transformation$^1$ (CT) OP means a symplectomorphism $f:M\to M$ on a symplectic manifold $(M,\omega)$ with Hamiltonian function $H:M\to\mathbb{R}$. Symplectomorphism are diffeomorphisms that respect the Poisson bracket $\{\cdot,\cdot\}$.
Darboux Theorem guarantees the existence of local Darboux/canonical coordinate charts. A symplectomorphism takes local Darboux/canonical coordinates into local Darboux/canonical coordinates.
Note that from a geometric point of view the important construct is the Poisson bracket itself. In contrast, the local Darboux/canonical coordinates are in principle dispensable.
Example: Let us here assume no explicit time dependence for simplicity. Then under a symplectomorphism $f$ a constant of motion $Q:M\to\mathbb{R}$ with $\{Q,H\} =0$ turn into a constant of motion $f^{\ast}Q:=Q\circ f$ with $\{f^{\ast}Q,f^{\ast}H\} =0$.
--
$^1$ For various definitions of a CT, see this Phys.SE post.
You're right. Sometimes it can be fine to take a non-canonical transformation and this is often done in practice (though not usually emphasized). I'll explain in an example from quantum mechanics. The analogy to classical mechanics can be seen by replacing the commutator with the Poisson bracket..
Consider a harmonic oscillator Hamiltonian
$$ H = \frac{1}{2} m \omega^2 X_0^2 + \frac{1}{2m}P_0^2 $$
We have that $X_0$ and $P_0$ satisfy the canonical commutation relations $[X_0,P_0] = i\hbar$. Suppose we want to think about the harmonic oscillator in phase space but with $X$ and $P$ on equal footing. We can rescale the variables in such a way that 1) the new variables have the same dimensions and 2) canonical commutation relations are preserved.
\begin{align} X' &= \sqrt{m\omega} X_0\\ P' &= \frac{1}{\sqrt{m\omega}} P_0 \end{align}
This gives us
\begin{align} H &= \frac{\omega}{2}(X'^2 + P'^2)\\ [X',P'] &= i\hbar \end{align}
This was a canonical transformation because it preserves the commutation relations.
Ok. Now I'm going to describe and perform a non-canonical transformation and afterwards I will motivate it more. Say that not only do we want $X$ and $P$ to have the same dimension but we actually want them to be dimensionless! In quantum mechanics there is $\hbar$ which conveniently has units of position times momentum (see canonical commutation relation). If we divide through by this we can make both dimensionless.
\begin{align} X &= \frac{1}{\sqrt{2\hbar}} X'\\ P &= \frac{1}{\sqrt{2\hbar}} P'\\ H &= \hbar \omega (X^2+P^2)\\ [X,P] &= \frac{i}{2} \end{align}
We see that this transformation is non-canonical because it changes the commutation relations. However, it is still useful because it brings the Hamiltonian into a very nice and recognizable form. First of all, we can see that
$$ X^2 + P^2 = a^{\dagger}a + \frac{1}{2} $$
Where $a$ and $a^{\dagger}$ are the annihilation and creation operators. The $\frac{1}{2}$ is unimportant for this discussion. A bit of work can show in this form we have
\begin{align} a &= X+iP\\ a^{\dagger} &= X-iP\\ X &= \frac{1}{2}(a^{\dagger}+a)\\ P &= \frac{i}{2}(a^{\dagger}-a) \end{align}
This transformation was really nice because now that $X$ and $P$ are unitless we can ALSO treat them on equal ground with the unitless bosonic creation and annihilation operators and it is very obvious that $X$ and $P$ are directly related to the real and imaginary parts (Hermitian and ant-Hermitian) parts of $a$.
All of this was to motivate that sometimes it is useful to do non-canonical transformations. Now a note as to why we might want to avoid non-canonical commutation relations. The answer is just that it introduces the possibility of making a mistake. Say we are looking for the equations of motion. We know the Heisenberg equation of motion for an operator is
$$ \dot{O} = -\frac{i}{\hbar}[O,H] $$
We can calculate the time evolution for $X$ as
\begin{align} \dot{X} &= -\frac{i}{\hbar}\hbar \omega [X,P^2] = -\frac{i}{\hbar} 2\frac{i}{2}\hbar \omega P = \omega P \end{align}
However, if we had given the Hamiltonian to our friend but forgotten to tell them that we are using $X$ and $P$ with $[X,P] = \frac{i}{2}$ then they may have accidentally used the canonical commutation relations and calculated
\begin{align} \dot{X} = -\frac{i}{\hbar}\hbar \omega [X,P^2] \rightarrow -\frac{i}{\hbar} \hbar \omega 2 i \hbar P= 2\hbar \omega P \end{align}
This would be incorrect. So basically, the canonical commutation relations are sort of like conventional commutation relations. You can use non-canonical transformations, and sometimes this may be useful, but allowing yourself to use non-canonical transformations means you have to do a better job keeping track of what variables you are using. In quantum mechanics this is something which is natural to do anyways so it is not really surprising or commented upon when such non-canonical variables are used. In classical mechanics it may be less useful and more strongly enforced that one should not use non-canonical coordinates so it is done less often.
See the Wikipedia on Optical Phase Space for an example where these sorts of $X$ and $P$ quadrature operators are used. Otherwise see any intro quantum optics book.
You are quite right, there is no objective reason why we should restrict ourselves to canonical transformations. As a matter of fact, the Hamilton equations can be formulated in a coordinate-independent way (by coordinates I mean the coordinates on phase space, i.e. the configuration variables and conjugate momenta), so that the very need of choosing coordinates does not arise. The coordinate-independent formulation is given by the following equation:
$$ \omega(X,\cdot)=dH $$
Here $H$ is your Hamiltonian and $dH$ is differential, $X$ is the tangent vector to the dynamical trajectory on phase space and $\omega$ is called a symplectic form. It is a 2-form (i.e. it is an antisymmetric tensor that accepts two vectors as its arguments) which moreover is closed, in the sense that $d\omega=0$, where the (exterior) differential of $\omega$ is defined in arbitrary coordinates $x^{I}$ (both configuration and momentum variables) as
$$ (d\omega)_{IJK}=\frac{\partial \omega_{JK}}{\partial x^{I}}+\frac{\partial \omega_{KI}}{\partial x^{J}}+\frac{\partial \omega_{IJ}}{\partial x^{K}} $$
The vector $X$ is defined in the same coordinates as
$$ X^{I}=\frac{dx^{I}}{dt} $$
Now, there is a theorem due to Darboux that says that if the phase space that you are considering admits a symplectic form, then locally one can find (non-unique) coordinates $(q^{i},p_{i})$ (the number of $I$ indices is double the number of $i$ indices) such that
$$ \omega=dq^{i}\wedge d{p}_{i} $$
where $\wedge$ is the wedge product, i.e. the antisymmetric part of the regular tensor product multiplied by two. Equivalently,
$$ \omega(X,Y)=X^{i}Y_{i}-X_{i}Y^{i} $$
where $X^{i}$ and $X_{i}$ are the components of the vector $X$ with respect to the given coordinates (the same goes for $Y$). If $X$ is the tangent vector to a curve on phase space, then
$$ X^{i}=\frac{dq^{i}}{dt}\qquad X_{i}=\frac{dp_{i}}{dt} $$
Canonical transformations are precisely those that keep $\omega$ invariant in the form given above.
There is also a theorem (rather, a definition) that endows the phase spaces that you derive from Lagrangian mechanics with a canonical symplectic form, so that the Darboux theorem always holds in such phase spaces (they are called cotangent bundles on configuration spaces).
Let's now see what the Hamilton equations given in the coordinate-free form translate into with respect to the Darboux coordinates. We have
$$ dH=\frac{\partial H}{\partial q^{i}}\,dq^{i}+\frac{\partial H}{\partial p_{i}}\,dp_{i} $$
and
$$ \omega(X,\cdot)=dq^{i}(X)dp_{i}-dp_{i}(X)dq^{i}=X^{i}dp_{i}-X_{i}dq^{i}=\frac{dq^{i}}{dt}\,dp_{i}-\frac{dp_{i}}{dt}\,dq^{i} $$
Hence
$$ \frac{dq^{i}}{dt}\,dp_{i}-\frac{dp_{i}}{dt}\,dq^{i}=\frac{\partial H}{\partial q^{i}}\,dq^{i}+\frac{\partial H}{\partial p_{i}}\,dp_{i} $$
or
$$ \frac{dq^{i}}{dt}=\frac{\partial H}{\partial p_{i}}\qquad \frac{dp_{i}}{dt}=-\frac{\partial H}{\partial q^{i}} $$
which are the usual Hamilton equations. However, since $\omega(X,\cdot)=dH$ does not make reference to any specific coordinates, the Darboux coordinates are as good as any other set of coordinates. Alternatively, one may just start from the usual form of the Hamilton equations and define arbitrary changes of coordinates, as you have figured out. So what is the use of coordinate transformations?
There are two answers to this question. The first one is that the Darboux coordinates, which are transformed into one another by canonical transformations, are the ones in which the Hamilton equations take arguably the simplest form. But you may not care about this simplification. Then there is a second, more relevant answer: it can be shown that any transformation on the configuration space (the $q^{i}$'s) which is a symmetry of the Lagrangian action of your system lifts automatically to a canonical transformation on the phase space of the Hamiltonian formulation, meaning that to any symmetry transformation of the configuration space is associated a canonical transformation of the phase space. Therefore, if you focus on symmetries rather than on the Hamilton equations, you find that the canonical transformations on phase space are the ones that realize a symmetry transformation on the system: solutions of the dynamics which are related to one another through a symmetry transformation are also related by a canonical transformation.
It's for the same reason that Einstein sought a covariant form for the equations of gravity. The equations should not be dependent on our coordinatisation of space. This he called general covariance.
In classical mechanics, canonical transformations leave invariant the Hamiltonian and this means it leaves invariant the time evoution of the system, and hence the physics f the system is still the same.
