If we have a system of PDE of the form:

$$\begin{cases} \dfrac{\partial y}{\partial t}(t,x)=D\Delta y+F(t,x,f(x),y) ,\ (t,x)\in (0,T)\times\Omega\\ \dfrac{\partial y}{\partial \nu}(t,x)=0,\ (t,x)\in (0,T)\times\partial\Omega \\ y(0,x)=y_{0}(x),\ x\in\Omega\end{cases}$$

with a unique solution $y=y^f$(mild solution). All derivatives are considered in the distributional sense (weak derivatives).Here $\Omega\subset\mathbb{R}^N$ is an open and bounded set with $C^1$ boundary, $y=(y_1,y_2,...,y_m)$, $F=(F_1,F_2,\dots,F_m)$, $D=diag(d_1,d_2,\dots, d_m),\ d_i\in\mathbb{R},\ \forall\ i=1,2,...,m$ and $f:\Omega\to\mathbb{R}$ a function in $L^1(\Omega)$. We consider the following approximation of the system:

$$\begin{cases} \dfrac{\partial y}{\partial t}(t,x)=D\Delta y+F(t,x,f_n(x),y) ,\ (t,x)\in (0,T)\times\Omega\\ \dfrac{\partial y}{\partial \nu}(t,x)=0,\ (t,x)\in (0,T)\times\partial\Omega \\ y(0,x)=y_{0}(x),\ x\in\Omega\end{cases}$$,

with the unique solution $y^{f_n}$.

How can we prove that if $f_n\to f$ in $L^1(\Omega)$, then $y_i^{f_n}\to y_i^f$ in $L^p((0,T)\times\Omega),\ \forall\ i=1,2,\dots,n$, for some $p\in [1,\infty)$? Are there some results in that sense? Maybe with a more general setting...?

Are there some estimates?

$$||y^{f_n}-y^f||\leq C||f_n-f||_{L^1(\Omega)}$$.