Lovelock and Horndeski theories of gravity are diffeomorphism-invariant theories with second-order equations of motion. A subset of these theories can be motivated by effective field theory considerations and hence they could describe strong-field deviations from general relativity. In particular, the effects of some Horndeski theories might be observable by present and future gravitational wave detectors. To study the dynamics of the theories using numerical simulations, they must satisfy some mathematical consistency properties. In this thesis, we establish two such properties for Lovelock and Horndeski theories.

In the first part of the thesis, we study the Cauchy problem for Lovelock and Horndeski theories. To demonstrate that the Cauchy problem for a theory of gravity is locally well-posed, it is sufficient to show that the gauge-fixed equations of motion are strongly hyperbolic. First, we use some numerical-relativity-inspired gauge conditions to write the equations of motion of weakly coupled cubic Horndeski theories in a strongly hyperbolic form. Next, this result is strengthened by proving that any weakly coupled Lovelock and Horndeski theory possesses a strongly hyperbolic formulation. This is achieved by introducing a novel class of "modified harmonic" gauge conditions and gauge-fixing procedures.

Another essential requirement on a theory of gravity is the possibility to choose initial data that represents astrophysically realistic systems. Some of the physically most interesting systems are approximately isolated systems and can be modelled by asymptotically flat spacetimes. The second part of the thesis discusses three methods to construct such initial data for a class of Horndeski theories. These methods are based on standard conformal techniques used in general relativity to write the constraint equations as a system of elliptic partial differential equations. It is shown that for a class of weakly coupled Horndeski theories, the conformally formulated constraint equations admit a well-posed boundary value problem on asymptotically Euclidean initial data surfaces.