Three Coupled PDE. Is my approach in the right way?

Question

I have the following three 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}

From the first and second equation i expressed $\frac{\partial \theta_h}{\partial x}$ and $\frac{\partial \theta_c}{\partial y}$ in terms of $\theta_w$. Then i substituted these in the third equation to yield

$$\lambda_h \frac{\partial^2 \theta_w}{\partial x^2} + \lambda_c V \frac{\partial^2 \theta_w}{\partial y^2} - (\beta_h+V\beta_c)\theta_w+(\beta_h\theta_h+V\beta_c\theta_c) = 0$$

The third PDE turns out to be a second order linear Elliptic PDE as $\lambda_h$,$\lambda_c$ and $V$ are all positive constants. I have reached a canonical form for this second order PDE. This PDE is defined on a rectangle with Neumann conditions. I plan to do the following next;

1. Calculate $\theta_w(x,y)$ from the second order PDE.
2. Plug them in the first two to obtain $\theta_h$ and $\theta_c$

Am i following a correct approach or is there any subtlety i am over-looking ?

Attempt

The boundary conditions for the problem are as follows:

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$$

After the suggestions from @Christoph here i have the following two linear third order differential equations:

\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}

Both these ODEs now need to be converted to individual Boundary value problems using the BC(s). On substituting

$$\theta_w(x,y) = e^{-\beta_h x} F(x) e^{-\beta_c y} G(y)$$ into the give BC(s), i arrive at the following

$$e^{-\beta_cy}F(0)G(y)=1$$ $$e^{-\beta_hx}F(x)G(0)=0$$ $$e^{-\beta_cy}G(y)[F'(0)-\beta_hF(0)]=0$$ $$e^{-\beta_cy}e^{-\beta_h}G(y)[F'(1)-\beta_hF(1)]=0$$ $$e^{-\beta_hx}F(x)[G'(0)-\beta_cG(0)]=0$$ $$e^{-\beta_hx}e^{-\beta_c}F(x)[G'(1)-\beta_cG(1)]=0$$

Following this (keeping in mind that exponential cannot attain a 0 value), i arrive at the following simplifications:

$$G(0)=0$$ $$G'(0)=0$$ $$\frac{G'(1)}{G(1)}=\beta_c$$ $$\frac{F'(0)}{F(0)}=\beta_h$$ $$\frac{F'(1)}{F(1)}=\beta_h$$

The ODEs are of third order and although i have six BC(s), on decoupling the BC(s) i get just 5. Am i misunderstanding something or is there some other way ?

Attempt 2

As @Christoph advised I made the following changes: $$\bar{{\theta_h}}(x,y):=\theta_h(x,y)-1$$ and the ansatz $$\theta_w(x,y)=e^{-\beta_hx}f(x)e^{-\beta_cy}g(y)$$ such that $F(x) := \int f(x) \, \mathrm{d}x$ and $G(y) := \int g(y) \, \mathrm{d}y$

The third order linear DEs we arrive at still remain the same. For figuring out the b.c.(s), the ansatz became: $$\theta_w(x,y)=e^{-\beta_hx}F'(x)e^{-\beta_cy}G'(y)$$

But the boundary conditions now take the following form

For $F$: $$F(0)=0$$ $$\frac{F''(0)}{F'(0)}=\beta_h$$ $$\frac{F''(1)}{F'(1)}=\beta_h$$

For $G$: $$G(0)=0$$ $$\frac{G''(0)}{G'(0)}=\beta_c$$ $$\frac{G''(1)}{G'(1)}=\beta_c$$

Now i have three b.c. (s) for each boundary value problem viz. $F$ and $G$. Each BVP (one each of $F$ and $G$) now involve one Dirichlet and two Robin type b.c.

Answer 1

Here are a few hints:

1. Solve the two first-order PDEs for $\theta_h, \theta_c$ as functions of $\theta_w$: \begin{eqnarray} \theta_h(x,y) &=& \beta_h e^{-\beta_h x} \int e^{\beta_h x} \theta_w(x,y) \, \mathrm{d}x,\\ \theta_c(x,y) &=& \beta_c e^{-\beta_c y} \int e^{\beta_c y} \theta_w(x,y) \, \mathrm{d}y. \end{eqnarray}

2. Eliminate $\theta_h, \theta_c$ in the second-order PDE to obtain the following equation for $\theta_w$: \begin{eqnarray} 0 &=& e^{-\beta_h x} \left( \lambda_h e^{\beta_h x} \frac{\partial^2 \theta_w}{\partial x^2} - \beta_h e^{\beta_h x} \theta_w + \beta_h^2 \int e^{\beta_h x} \theta_w \, \mathrm{d}x \right) +\\ && + V e^{-\beta_c y} \left( \lambda_c e^{\beta_c y} \frac{\partial^2 \theta_w}{\partial y^2} - \beta_c e^{\beta_c y} \theta_w + \beta_c^2 \int e^{\beta_c y} \theta_w \, \mathrm{d}y \right). \end{eqnarray}

3. Use separation of variables with the ansatz $\theta_w(x,y) = e^{-\beta_h x} f(x) e^{-\beta_c y} g(y)$. You should obtain two linear third-order ODEs with constant coefficients for $F(x) := \int f(x) \, \mathrm{d}x$ and $G(y) := \int g(y) \, \mathrm{d}y$: \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} with some separation constant $\mu \in \mathbb{R}$.

4. From the boundary conditions on $\theta_h, \theta_c, \theta_w$ you should obtain conditions on $F$, $G$. Once the BVPs with the third-order ODEs are solved, compute $f \equiv F'$, $g \equiv G'$.