We introduce a novel constructive approach to define time evolution of critical points of an energy functional. Our procedure, which is different from other more established approaches based on viscosity approximations in infinite dimension, is prone to efficient and consistent numerical implementations, and allows for an existence proof under very general assumptions. We consider in particular rather nonsmooth and nonconvex energy functionals, provided the domain of the energy is finite dimensional. Nevertheless, in the infinite dimensional case study of a cohesive fracture model, we prove a consistency theorem of a discrete-to-continuum limit. We show that a quasistatic evolution can be indeed recovered as a limit of evolutions of critical points of finite dimensional discretizations of the energy, constructed according to our scheme. To illustrate the results, we provide several numerical experiments both in one and two dimensions. These agree with the crack initiation criterion, which states that a fracture appears only when the stress overcomes a certain threshold, depending on the material.Moreover, we say that a measurable function v : [0, T ] → E is an approximable quasistatic evolution with initial condition v 0 and constraint f , if for every t ∈ [0, T ] there exists a sequence δ k → 0 + and a sequence (v δ k ) k∈N of discrete quasistatic evolutions with time step δ k , initial condition v 0 , and constraint f , such thatWe are now ready to state the main results of the paper (see Theorem 2.15 and Theorem 2.16).Dissipative systems. (0)). Under suitable assumptions on J, ψ, A, and f (see Theorem 2.15) we prove that thereis an approximable quasistatic evolution with initial condition v 0 and constraint f ;(C) Energy inequality: the function s → q(s),ḟ (s) F belongs to L 1 (0, T ) andis defined by (2.9) and ·, · F denotes the duality product in F.Any function v ∈ BV ([0, T ]; E) satisfying (A), (B), and (C) is a rate independent evolution. We also remark that, as a consequence of the local stability (B), the energy inequality becomes an equality in all the nontrivial intervals where the solution v(t) is absolutely continuous (see Theorem 2.18, where also the nondissipative case is treated). In addition, in such intervals the doubly nonlinear inclusion A * q(t) ∈ ∂J(v(t)) + ∂ψ(v(t)) is satisfied. It is however a well-known fact that, due to nonconvexity, our solutions can in general develop time discontinuities, where additional dissipation appears (see [27]). Non dissipative systems. When there is no dissipation we can still prove an existence result, although the evolution obtained is in general expected to be less regular in time. This is due to the fact that the absence of dissipation causes loss of compactness, since simple estimates of the total variation of the approximating solutions are now missing. This is undoubtedly a point of great interest in the analysis of the degenerate case (see also [2] for an abstract approach in the unconstrained setting). To compensate the loss of compactness, we need ...