In this paper, we discuss the learning and discovery problem for the dynamical systems described through stable evolutionary Partial Differential Equations (PDEs). The main idea is to employ a suitable learning approach for creating a map from boundary conditions to the corresponding output. More precisely, in order to accurately uncover the evolutionary PDE dynamics, we propose a scheme that employs large-scale system identification to construct such a map using sufficiently informative measurements. Accordingly, we first develop a scalable implementation for the subspace identification method, enforcing stability on the identified system. To this end, numerical optimization techniques such as coordinate descent, randomized singular value decomposition, and large-scale semidefinite programming are employed. The performance and complexity of the resulting scheme are discussed and demonstrated through numerical experiments on generic identification examples. Following this, we validate the effectiveness of the proposed approach on an example of a stable evolutionary partial differential equation. The numerical results confirm the efficacy of the proposed learning scheme.