In this work, we propose to learn time integration schemes based on neural networks which satisfy three distinct sets of mathematical constraints, i.e., unconstrained, 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 prediction 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 prediction error on coarse grids, with an overall better performance for the semi-constrained model. On a grid $4×$ 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 substantial improvement in phase prediction for the wave equation. On a $32×$ 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.