Creation of the weak formulation is a standard exercise. We assume $\Omega \subset \mathbb R^n$.
Indeed, suppose that $$
\Delta u(x) - v(x) = f(x) \quad ,\quad \forall \ x \in \Omega \tag{1}
$$
Let $\phi$ be a function in $\Omega$ whose regularity will be decided based on desired manipulations. Then, assuming that $(1)$ holds, we multiply by $\phi$ (note : while writing integrals, I omit the variable $x$ and the notation $dx$ to focus on the integrand. This is common practice in PDE):
$$
\phi\Delta u - \phi v = \phi f \implies \int_{\Omega} \phi \Delta u - \int_{\Omega} \phi v
= \int_{\Omega} \phi f$$
Let's deal with just the first term : $\int_{\Omega} \phi \Delta u$, the rest will remain as they are. We write by extension by zero (allow $\phi$ to be zero on $\partial \Omega$, so the function remains integrable) :
$$
\int_{\Omega} \phi \Delta u = \int_{\mathbb R^n} \phi \Delta u
$$
Under suitable conditions, Fubini's theorem and integration-by-parts apply and we can write :
$$
\int_{\mathbb R^n} \phi \Delta u = \int_{\mathbb R^{n-1}} \left[\int_{\mathbb R} \phi \Delta u \right] = -\int_{\mathbb R^{n-1}} \int_{\mathbb R} \nabla \phi \nabla u = -\int_{\mathbb R^n} \nabla \phi \nabla u = -\int_{\Omega} \nabla \phi \nabla u
$$
where the second equality follows from integration-by-parts and the fact that $\phi$ is compactly supported hence the derivatives of it vanish at infinity.
Which means we land up with the weak formulation for this particular equation :
$$
-\int_{\Omega} \nabla \phi \nabla u - \int_{\Omega} \phi v = \int_{\Omega} \phi f
$$
The conditions are quite self-evident : we require (weakly) $\nabla u , \nabla \phi$ to exist and be integrable , and furthermore vanish on the boundary. It is only natural to then call for $\phi,u \in H_0^1(\Omega)$.
The other equation is solved similarly, and results in
$$
\int_{\Omega} \nabla \phi \nabla v + \int_{\Omega} \phi u = \int_{\Omega} \phi g
$$
for all $\phi \in H_0^1(\Omega)$. Note that $v \in H_0^1(\Omega)$ from here as well.
To create our bilinear form is now extremely clear : we must "group" the two equations together. Define :
$$
B((\phi,\phi'),(u,v)) = \int_{\Omega} -\nabla \phi \nabla u - \int_{\Omega} \phi v + \int_{\Omega} \nabla \phi' \nabla v + \int_{\Omega} \phi' u
$$
so $B$ is a bilinear form on $H_0^1(\Omega) \times H_0^1(\Omega)$. Define :
$$
J : H_0^1(\Omega) \times H_0^1(\Omega) \to \mathbb R \quad ; \quad J(\phi,\phi') = \int_{\Omega} \phi f + \int_{\Omega} \phi' g
$$
this is a well defined linear functional. Our equation then becomes : find $(u,v)$ such that $$
B((\phi,\phi'),(u,v)) = J(\phi,\phi')
$$
for all $(\phi,\phi') \in H_0^1(\Omega)\times H_0^1(\Omega)$. We are firmly in Lax-Milgram territory ; the LHS is a bilinear form in $H_0^1(\Omega) \times H_0^1(\Omega)$ and the RHS is a linear functional on the same space.
To apply the Lax-Milgram theorem, whose direct consequence is the existence and uniqueness of $u,v$ satisfying the required equation (as desired), we must prove that $B$ is continuous and coercive, and $J$ is continuous.
Continuity of either one is fairly obvious :
$$
|B((\phi,\phi'),(u,v))| = \left|\int_{\Omega} -\nabla \phi \nabla u - \int_{\Omega} \phi v + \int_{\Omega} \nabla \phi' \nabla v + \int_{\Omega} \phi' u \right| \\
\leq \int_{\Omega} |\nabla \phi| |\nabla u| + \int_{\Omega} |\phi| |v| + \int_{\Omega} |\nabla \phi'| |\nabla v| + \int_{\Omega} |\phi'||u| \\
\leq \|\phi\|\|u\| + \|\phi\|\|v\| + \|\phi'\|\|u\| + \|\phi'\|\|v\| \\
\leq (\|(\phi,\phi')\|_{H_0^1 \times H_0^1}) (\|(u,v)\|_{H_0^1 \times H_0^1})
$$
where on the third line all norms are in $H_0^1$. I leave you to figure out the case for $J$.
The coercivity of $B$ follows from Poincare's inequality. Indeed, by Poincare's inequality, which can be used since $\Omega$ is bounded, we know that there exists a positive constant $C>0$ such that $\int_{\Omega} \nabla h \nabla h \geq C \|h\|^2_{H_0^1}$ for every $h \in H_0^1(\Omega)$.
This helps massively, since :
$$
B((u,v),(u,v)) = \int_{\Omega} \nabla u \nabla u + \int_{\Omega} v \cdot v + \int_{\Omega} \nabla v \nabla v + \int_{\Omega} u \cdot u \geq C(\|u\|_{H_0^1} + \|v\|_{H_0^1}) = C (\|(u,v)\|_{H_0^1 \times H_0^1})
$$
because the terms $u \cdot u , v \cdot v$ are non-negative, and for the rest the Poincare inequality apply (note : I'm putting the dots to indicate that it is a dot product : it was one earlier as well between the $\phi$s and $u,v$ but because they were different variables I did not use the dot)
Thus , we may use Lax-Milgram and complete the problem.
