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?
By using the answer here with $m=\|\varphi\|_\infty$, we can assume that there exist a sequence $\varphi_k\in C_0^\infty$ such that $\varphi_k\to \varphi$ in $W^{1,p}$ and $\|\varphi_k\|_\infty\leq M$ where $M$ is a positive constant.
We use Theorem 4.9. from Brezis to conclude without loss of generality that $\varphi_k\to\varphi$ almost everywhere, hence $a_0(x,u,Du)\varphi_k\to a_0(x,u,Du)\varphi$ almost everywhere and $|a_0(x,u,Du)\varphi_k|\leq M|a_0(x,u,Du)|\in L^1$, therefore we can apply Lebesgue theorem to conclude that
$$\tag{1}\int_\Omega a_0(x,u,Du)\varphi_k\to\int_\Omega a_0(x,u,Du)\varphi$$
On the other hand, the functions $a_i(x,u,Du)$ defines bounded linear functionals in $W_0^{1,p}$, which implies that $$\tag{2}\int_\Omega a(x,u,Du)D\varphi_k\to\int_\Omega a(x,u,Du)D\varphi$$
To conclude, we combine $(1)$ and $(2)$ to get the desired equality.
Remark: Just for the sake of clarity, I will prove here that the operator defined by @user98130 in his answer is continuous from $W_0^{1,p}$ to $W_0^{1,p}$.
First note that $\psi$ is Lipschitz, therefore, we can assume that $$|\psi(x)-\psi(y)|\leq C |x-y|,\ \forall\ x,y\in\mathbb{R}$$
Now, suppose that $u_k\to u$ in $W_0^{1,p}$. Note that $$\int_\Omega |\psi(u_k(x))-\psi(u(x)|^p\leq\int_\Omega C^p |u_k(x)-u(x)|^p\tag{3}$$
On the other hand $$\left\|\psi'(u_k)\frac{\partial u_k}{\partial x_i}-\psi'(u)\frac{\partial u}{\partial x_i}\right\|_p\leq \left\|(\psi'(u_k)-\psi'(u))\frac{\partial u_k}{\partial x_i}\right\|_p+\left\|\psi'(u)\left(\frac{\partial u_k}{\partial x_i}-\frac{\partial u}{\partial x_i}\right)\right\|_p\tag{4}$$
We apply Theorem 4.9. of Brezis and then Lebesgue Theorem to conclude that the first term from the right hand side of $(4)$ converges to $0$. To this end, note that by Theorem 4.9. from Brezis's book, we can assume that $u_k\to u$ almost everywhere, which implies that $\psi'(u_k)\to \psi'(u)$ almost everywhere. Moreover $\psi'(u_k)\leq M$, hence, by Lebesgue Theorem, $\psi'(u_k)\to \psi'(u)$ in $L^q$ for all $q\in [1,\infty)$. To conclude, just use Holder inquality.
The second term converges to $0$, because $\psi'$ is bounded and $u_k\to u$ in $W^{1,p}$.
We conclude from $(3)$ and $(4)$ the continuity of the composition operator. Another way to prove the continuity is: Let $T:W_0^{1,p}\to W_0^{1,p}$ be the composition operator, i.e. $T(u)=\psi\circ u$. You can verify that $T$ is Frechet differentiable and that $$T'(u)v=\psi'(u)v$$
Moreover, although this is not necessary, you can also verify that $T$ is $C^1$..
Remark 1: You can also use this approach to prove that $\varphi_k$ can be choose to be bounded.