  In this paper, we study flexibility of the weak solutions to the
Monge-Amp\`ere system (MA) via convex integration. This system of Pdes is an
extension of the Monge-Amp\`ere equation in $d=2$ dimensions, naturally arising
from the prescribed curvature problem, and closely related to the classical
problem of isometric immersions (II).
  Our main result achieves density, in the set of subsolutions, of the H\"older
regular $\mathcal{C}^{1,\alpha}$ solutions to the weak formulation (VK) of
(MA), for all $\alpha<\frac{1}{1+d(d+1)/k}$ where $d$ is an arbitrary dimension
and $k$ is an arbitrary codimension of the problem. This result seems to be
optimal, from the technical viewpoint, for the corrugation-based convex
integration scheme. In particular, it covers the codimension interval $k\in
\big(1, d(d+1)\big)$, so far uncharted even for the system (II), since the
regularity $\mathcal{C}^{1,\alpha}$ with any $\alpha <1$ proved by K\"allen,
strictly requires a large codimension $k\geq d(d+1)$. We also reproduce
K\"allen's result in the context of (MA). At $k=1$, our result agrees with the
regularity $\mathcal{C}^{1,\alpha}$ for (II) with any $\alpha
<\frac{1}{1+d(d+1)}$, proved by Conti, Delellis and Szekelyhidi. Finally, our
results extend the initial findings by the author and Pakzad for (MA) and $d=2,
k=1$.
  As an application of our results for (VK), we derive an energy scaling bound
in the quantitative immersability of Riemannian metrics, for nonlinear energy
functionals modelled on the energies of deformations of thin prestrained films
in the nonlinear elasticity.
