A complicated boundary value problem leading to third order Eigen sytem [Help in continuing forward] [EDITED]

Question

I have the following elliptic PDE (describing temperature in a plate, w in thermal contact with two fluids h and c): $$\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 \tag 1$$

pde = Subscript[λ, h] D[θ[x, y], x, x] + Subscript[λ, c] V D[θ[x, y], y, y] - (Subscript[β, h] + V Subscript[β, c] [Theta][x, y] + (Subscript[β, h] θh[x, y] + V Subscript[β, c] θc[x, y]) == 0

with the boundary conditions ($\beta_h, \beta_c, V, \lambda_h, \lambda_c$ are constants)

$$\theta_w(0,y)=1, \theta_w(x,0)=0$$ $$\frac{\partial \theta_w(1,y)}{\partial x}=\frac{\partial \theta_w(x,1)}{\partial y}=0$$

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

It is also known that $\theta_h(0,y)=1$ and $\theta_c(x,0)=0$.

Trying with separation of variables:

Substituting $\theta_h$ and $\theta_c$ in $(1)$ and using the ansatz $\theta_w(x,y) = e^{-\beta_h x} f(x) e^{-\beta_c y} g(y)$ we 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, \tag 2\\ 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,\tag 3 \end{eqnarray} with some separation constant $\mu \in \mathbb{R}$.

The boundary conditions take the form:

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

In order to get another homogeneous condition of $F$, we use $\theta_{h1}(x,y):=\theta_h (x,y)-1$. This does not change Eq. $(2)$ or $(3)$ but gives us another condition on $F$ i.e. $F'(0)=0$

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

Can anyone suggest a way forward in this problem ? I know I need to find the Eigenvalues $\mu$ and use the non-homogeneous b.c. and orthogonality to complete the solution, Krazug has already provided a method to find the $\mu$ values in an answer below. I need to understand how to use these with the $F$ equation to complete the $\theta_w$ solution

---

CONTEXT The integro-differential equation from which the above problem arises is $D$:

I had a system of three PDEs $$\frac{\partial \theta_h}{\partial x}+\beta_h (\theta_h-\theta_w) = 0 \tag A$$

$$\frac{\partial \theta_c}{\partial y} + \beta_c (\theta_c-\theta_w) = 0 \tag B$$

$$\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 \tag C$$

On eliminating $\theta_h$ and $\theta_c$ from $C$ I reach $$\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^2 e^{-\beta_h x} \int e^{\beta_h x} \theta_w(x,y) \mathrm{d}x + \beta_c^2 e^{-\beta_c y}\int e^{\beta_c y} \theta_w(x,y)\mathrm{d}y = 0 \tag D$$

The BCs are

$$\theta_w(0,y)=1, \theta_w(x,0)=0$$ $$\frac{\partial \theta_w(1,y)}{\partial x}=\frac{\partial \theta_w(x,1)}{\partial y}=0$$ and $\theta_h(0,y)=1$, $\theta_c(x,0)=0$.

Answer 1

This is not a full answer to the question!

I have a package which helps solve eigenvalue BVPs by calculating the Evans function, an analytic function whose roots correspond to the eigenvalues. Some details are available at these two questions, or this PDF. Or search for CompoundMatrixMethod to see my previous answers here. It comes much better than the determinant method for finding roots accurately.

Install the package (also available on my github page):

Needs["PacletManager`"]
    PacletInstall["CompoundMatrixMethod", 
    "Site" -> "http://raw.githubusercontent.com/paclets/Repository/master"]

Load the package and setup the system:

Needs["CompoundMatrixMethod`"]

eqn = {V λc G'''[x] - 2 V λc βc G''[x] + ((λc βc - 1) V βc + μ) G'[x] + V βc^2 G[x] == 0};
bcs = {G[0] == 0, G'[0] == 0, G''[1] == βc G'[1]};
sys = ToMatrixSystem[eqn, bcs, {G}, {x, 0, 1}, μ] /. βc ->
        0.921 /. λc -> 1.775*10^-4;

Now we can evaluate the Evans function at a given value of $\mu$ and $V$ with e.g. Evans[μ, sys /. V -> 1]. Roots of this function correspond to eigenvalues of that third order equation.

Plot[{Evans[μ, sys /. V -> 1], Evans[μ, sys /. V -> 2]}, {μ, 2, 3}]

You can see there is a periodic set of eigenvalues coming out for positive values of $\mu$, and it is worth noting that the eigenvalues for $V=2$ are exactly twice those for $V=1$.

FindRoot[Evans[μ, sys /. V -> 1], {μ, 2}]
{μ -> 1.97219}