  Within the landscape of modified theories of gravity, progress in
understanding the behaviour of, and developing tests for, screening mechanisms
has been hindered by the complexity of the field equations involved, which are
nonlinear in nature and characterised by a large hierarchy of scales. This is
especially true of Vainshtein screening, where the fifth force is suppressed by
high-order derivative terms which dominate within a radius much larger than the
size of the source, known as the Vainshtein radius.
  In this work, we present the numerical code $\varphi$enics, building on the
FEniCS library, to solve the full equations of motion from two theories of
interest for screening: a model containing high-order derivative operators in
the equation of motion and one characterised by nonlinear self-interactions in
two coupled scalar fields. We also include functionalities that allow the
computation of higher-order operators of the scalar fields in post-processing,
enabling us to check that the profiles we find are consistent solutions within
the effective field theory. These two examples illustrate the different
challenges experienced when trying to simulate such theories numerically, and
we show how these are addressed within this code. The examples in this paper
assume spherical symmetry, but the techniques may be straightforwardly
generalised to asymmetric configurations. This article therefore also provides
a worked example of how the finite element method can be employed to solve the
screened equations of motion. $\varphi$enics is publicly available and can be
adapted to solve other theories of screening.
