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$