Let $\Omega\subset \mathbb{R}^n, \ n\geq 2$ be a bounded Lipschitz domain. Consider the following inhomogenous system (in divergence form) subject to zero Dirichlet boundary conditions: \begin{equation} -\sum_{i=1}^{n}D_ia_i(x, u, Du)+a_0(x, u, Du)=0\quad\text{in } \Omega, \end{equation}where the vector field $a=(a_1, \dots, a_n)$ and the inhomogeneity $a_0$ satisfy the following (natural) growth condition: \begin{align} |a_i(x, u, Du)|&\leq K(1+|z|^{p-1}), \quad i=1, \dots, n,\\ |a_0(x, u, Du)|&\leq K_0(1+|z|^p) \end{align} for all $x\in \Omega,\ u\in \mathbb{R}^N$ and $z\in \mathbb{R}^{nN}$. We suppose that $a_i (i=0, 1, \dots, n)$ are Caratheodory functions.
A weak solution to the above system is a function $u\in W^{1,\ p}(\Omega, \mathbb{R}^N)$ satisfies \begin{equation*} \int_{\Omega}a(x, u, Du)\cdot D\varphi\ \mathrm{d}x+\int_{\Omega}a_0(x, u, Du)\cdot\varphi\ \mathrm{d}x=0\quad\forall\ \varphi\in C_c^{\infty}(\Omega, \mathbb{R}^N). \end{equation*}
In various papers/books about this topic, see for example Frehse and Beck (2013) "Regular and irregular solutions for a class of elliptic systems in critical dimensions", it is said that by approximation, the above identity remains true for all $\varphi\in W_0^{1,\ p}(\Omega, \mathbb{R}^N)\cap L^{\infty}(\Omega, \mathbb{R}^N)$.
Does anyone know of this approximation argument?