Let $x \in \Omega$ ($\Omega \subset \mathbb{R}^2$ is a bounded set with a sufficiently smooth boundary) and let us suppose that we have a PDE defined as
$$\begin{cases} \dfrac{\partial u(t,x)}{\partial t} &= A(u(t,x)) + F(u(t,x)),\\[10pt] u(t,x)|_{\partial \Omega} &= 0, \end{cases}$$
in which $A$ is a linear operator (which does not depend on either $t$ or $x$ so it is basically a matrix) and $F$ is a nonlinear operator with some nice properties (bounded, Lipschitz property, etc.) which describes the spatial dependence. We know that this equation has a unique solution.
What are the usual techniques to show that the solution of this problem depends continuously on the right-hand side, meaning that if we consider the modified problem
$$\begin{cases} \dfrac{\partial u_{\varepsilon}(t,x)}{\partial t} & = A(u_{\varepsilon}(t,x)) + F(u_{\varepsilon}t,x)) + \varepsilon,\\[10pt] u_{\varepsilon}(t,x)|_{\partial \Omega} &= 0, \end{cases}$$ then $\Vert u - u_{\varepsilon} \Vert \rightarrow 0$ as $\varepsilon \rightarrow 0$?
I have tried to use the variation of parameters formula (the operator semigroup version) but to no avail. I would be grateful even for some references since usually "continuous dependence" means the one considering the initial value, and not the right-hand side one.
Start by remarking that $u_\epsilon -u$ verifies $$ \frac{(\partial u_\epsilon(t,x)-u(t, x))}{\partial t}=A(u(t, x)-u_\epsilon(t,x))+F(u(t, x))+F(u_\epsilon(t, x))-F(u(t, x)))+\epsilon $$ Now multiplying the previous equation by $u_\epsilon -u$ and integrating over $\Omega$ we obtain \begin{align} \dfrac{1}{2} \dfrac{d}{dt} \int_\Omega \vert u_\epsilon(t,x)-u(t, x) \vert^2&=\int_\Omega A (u_\epsilon(t,x)-u(t, x))^2\\ &+\int_\Omega F(u_\epsilon(t, x))-F(u(t, x))) (u_\epsilon(t, x)-u(t,x))\\ &+ \int_\Omega \epsilon (u_\epsilon(t,x)-u(t, x)) \end{align} Now, you use the assumptions you have on $A$ (Matrix) and $F$ (L-Lipschitz), the Cauchy-Schwartz inequality as well as the inequality $ab \leq \dfrac{1}{2} (a^2+b^2) \; \forall a,b \in \mathbb{R}$ to eventually obtain \begin{align} \dfrac{1}{2} \dfrac{d}{dt} \int_\Omega \vert u_\epsilon(t,x)-u(t, x) \vert^2&\leq \Vert A \Vert \Vert u_\epsilon-u\Vert_{L^2(\Omega)}^2\\ &+L \Vert u_\epsilon-u \Vert_{L^2(\Omega)}^2\\ &+ \dfrac{\epsilon^2}{2} + \dfrac{1}{2} \Vert u_\epsilon -u \Vert^2_{L^2(\Omega)} \end{align} Now use Gronwall inequality to obtain $$ \Vert u_\epsilon-u \Vert^2_{L^2(\Omega)} \leq exp\big(2t \big( \Vert A \Vert + L +1 \big)\Bigg( \bigg(\Vert u_\epsilon(0,.)-u(0,.) \Vert^2_{L^2(\Omega)} \Big)+ \epsilon^2 t \Bigg) $$ I'm assuming that for the initial conditions you have $u_\epsilon(0,x)=u(0,x)+\epsilon$ as well, so by substitution you get $$ \Vert u_\epsilon-u \Vert^2_{L^2(\Omega)} \leq exp\big(2t \big( \Vert A \Vert + L +1 \big)\Bigg( \bigg(\epsilon^2 meas(\Omega) \Big)+ \epsilon^2 t \Bigg) $$ Finally you can deduce that $\underset{\epsilon \rightarrow 0}{\lim} \Vert u_\epsilon-u \Vert^2_{L^2(\Omega)}=0.$