The purpose of this paper is threefold: (i) to introduce and study the Euler-Lagrange prolongations of flatness PDEs solutions (best approximation of flatness) via associated least squares Lagrangian densities and integral functionals on Riemannian manifolds; (ii) to analyze some decomposable multivariate dynamics represented by Euler-Lagrange PDEs of least squares Lagrangians generated by flatness PDEs and Riemannian metrics; (iii) to give examples of explicit flat extremals and non-flat approximations.