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