We theoretically study the quench dynamics for an isolated Heisenberg spin chain with a random on-site magnetic field, which is one of the paradigmatic models of a many-body localization transition. We use the time-dependent variational principle as applied to matrix product states, which allows us to controllably study chains of a length up to L=100 spins, i.e., much larger than L≃20 that can be treated via exact diagonalization. For the analysis of the data, three complementary approaches are used: (i) determination of the exponent β which characterizes the power-law decay of the antiferromagnetic imbalance with time; (ii) similar determination of the exponent βΛ which characterizes the decay of a Schmidt gap in the entanglement spectrum; and (iii) machine learning with the use, as an input, of the time dependence of the spin densities in the whole chain. We find that the consideration of the larger system sizes substantially increases the estimate for the critical disorder Wc that separates the ergodic and many-body localized regimes, compared to the values of Wc in the literature. On the ergodic side of the transition, there is a broad interval of the strength of disorder with slow subdiffusive transport. In this regime, the exponents β and βΛ increase, with increasing L, for relatively small L but saturate for L≃50, indicating that these slow power laws survive in the thermodynamic limit. From a technical perspective, we develop an adaptation of the "learning by confusion" machine-learning approach that can determine Wc.