  We consider the squared singular values of the product of $M$ standard
complex Gaussian matrices. Since the squared singular values form a
determinantal point process with a particular Meijer G-function kernel, the gap
probabilities are given by a Fredholm determinant based on this kernel. It was
shown by Strahov \cite{St14} that a hard edge scaling limit of the gap
probabilities is described by Hamiltonian differential equations which can be
formulated as an isomonodromic deformation system similar to the theory of the
Kyoto school. We generalize this result to the case of finite matrices by first
finding a representation of the finite kernel in integrable form. As a result
we obtain the Hamiltonian structure for a finite size matrices and formulate it
in terms of a $(M+1) \times (M+1)$ matrix Schlesinger system. The case $M=1$
reproduces the Tracy and Widom theory which results in the Painlev\'e V
equation for the $(0,s)$ gap probability. Some integrals of motion for $M = 2$
are identified, and a coupled system of differential equations in two unknowns
is presented which uniquely determines the corresponding $(0,s)$ gap
probability.
