I am not specifically asking for a solution, but any reference on any method i could read about would be a big help. This I clarify as I am aware of the fact that MathOverflow only deals with research level questions. A starting point is all I aim to get from posting this question here
I have the following coupled PDEs: \begin{eqnarray} \frac{\partial \theta_h}{\partial x} + \beta_h (\theta_h - \theta_w) &=& 0,\\ \frac{\partial \theta_c}{\partial y} + \beta_c (\theta_c - \theta_w) &=& 0,\\ \lambda_h \frac{\partial^2 \theta_w}{\partial x^2} + \lambda_c V \frac{\partial^2 \theta_w}{\partial y^2} - \frac{\partial \theta_h}{\partial x} - V\frac{\partial \theta_c}{\partial y} &=& 0 \end{eqnarray} The boundary conditions are :
The PDE needs to be solved on a rectangular region where $x$ varies between $0$ to $1$ and $y$ varies between $0$ to $1$.
$$\frac{\partial \theta_w(0,y)}{\partial x}=\frac{\partial \theta_w(1,y)}{\partial x}=0 $$
$$\frac{\partial \theta_w(x,0)}{\partial y}=\frac{\partial \theta_w(x,1)}{\partial y}=0 $$
$$\theta_h(0,y)=1 $$$$\theta_c(x,0)=0$$
$\beta_h,\beta_c,\lambda_h,\lambda_c,V$ are all constants $>0$
Attempt
The third equation can be written as $$\lambda_h \frac{\partial^2 \theta_w}{\partial x^2} + \lambda_c V \frac{\partial^2 \theta_w}{\partial y^2} = \frac{\partial \theta_h}{\partial x} + V\frac{\partial \theta_c}{\partial y} $$
Then, from the first two equations the following can be written :-
$$\lambda_h \frac{\partial^2 \theta_w}{\partial x^2} + \lambda_c V \frac{\partial^2 \theta_w}{\partial y^2} =\beta_h e^{-\beta_h x} \int e^{\beta_h x} \theta_w(x,y) \, \mathrm{d}x + \beta_c e^{-\beta_c y} \int e^{\beta_c y} \theta_w(x,y) \, \mathrm{d}y$$
This resulting equation looks a lot like Laplace equation with integral source terms (if the LHS is taken into canonical form).
What kind of procedure should I take up to solve this problem? Is there any standard problem type that corresponds to this situation.
Attempt taken forward
It can be seen that the resulting equation can be variable separated if we consider the following ansatz:
$\theta_w(x,y) = e^{-\beta_h x} f(x) e^{-\beta_c y} g(y)$
with $F(x) := \int f(x) \, \mathrm{d}x$ and $G(y) := \int g(y) \, \mathrm{d}y$
to get:
\begin{eqnarray} \lambda_h F''' - 2 \lambda_h \beta_h F'' + \left( (\lambda_h \beta_h - 1) \beta_h - \mu \right) F' + \beta_h^2 F &=& 0,\\ V \lambda_c G''' - 2 V \lambda_c \beta_c G'' + \left( (\lambda_c \beta_c - 1) V \beta_c + \mu \right) G' + V \beta_c^2 G &=& 0, \end{eqnarray}
for some separation constant $\mu$