In this work, we propose to learn time integration schemes based on neu-ral networks which satisfy three distinct sets of mathematical constraints, i.e., uncon-strained, semi-constrained with the root condition, and fully-constrained with both root and consistency conditions. We focus on the learning of 3-step linear multistep methods, which we subsequently applied to solve three model PDEs, i.e., the heat equation, the wave equation, and the Burgers' equation. The results show that the pre-diction error of the learned fully-constrained scheme is close to that of the Runge-Kutta method and Adams-Bashforth method. Compared to the traditional methods, the learned unconstrained and semi-constrained schemes significantly reduce the predic-tion error on coarse grids, with an overall better performance for the semi-constrained model. On a grid 4x coarser than the reference, the mean square error (MSE) shows a reduction of up to two order of magnitude for some heat equation cases, and a sub-stantial improvement in phase prediction for the wave equation. Onta 32x coarser grid, the MSE for the Burgers' equation can be reduced by up to 40% to 45%. Tests on the two-dimensional heat equation and wave equation show a good generalization capability for the constant optimized coefficients.