Can this boundary value problem leading to a partio-integral DE be solved using finite Fourier or any integral transform?

I asked the three-dimensional version of this problem here which lead to a trivial solution. I have now tried it in 2-D

$$frac{partial theta_h}{partial x} + b_h (theta_h – theta_w) = 0, tag 1\$$
$$frac{partial theta_c}{partial y} + b_c (theta_c – theta_w) = 0,tag 2\$$
$$lambda_h frac{partial^2 theta_w}{partial x^2} + lambda_c V frac{partial^2 theta_w}{partial y^2} + b_h (theta_h – theta_w) + V b_c (theta_c – theta_w) = 0 tag 3$$

It is known that $theta_h(0,y)=theta_{hi}$ and $theta_c(x,0)=theta_{ci}$. We know that $b_h,b_c,V,theta_{hi}, theta_{ci}$ are constants greater than 0. Using this information and some manipulation on $(1),(2)$ we can write $(3)$ as follows:

$$lambda_h frac{partial^2 theta_w}{partial x^2} + lambda_c V frac{partial^2 theta_w}{partial y^2} +( -b_h – V b_c )theta_w +b_h^2 e^{-b_h x} int e^{b_h x} theta_w(x,y) mathrm{d}x + Vb_c^2 e^{-b_c y}int e^{b_c y} theta_w(x,y)mathrm{d}y = 0$$

which can be manipulated to:

$$lambda_h frac{partial^2 theta_w}{partial x^2} + lambda_c V frac{partial^2 theta_w}{partial y^2} +( -b_h – V b_c )theta_w + b_he^{-b_h x}bigg(theta_{hi}+b_hint_0^x e^{b_h s}theta_w (s,y)mathrm{d}sbigg) + Vb_ce^{-b_c y}bigg(theta_{ci}+b_cint_0^y e^{b_c s}theta_w (x,s)mathrm{d}sbigg)=0 tag 4$$

$(4)$ is subjected to the following boundary conditions:

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

In mathematica code $(4)$ and the boundary conditions are:

eq = λh D[θw[x, y], {x, 2}] + λc V D[θw[x, y], {y, 2}] + (-bh - V bc) θw[x, y] + bh E^(-bh x) (θhi + bh Integrate[E^(bh s) θw[s, y], {s, 0, x} ]) + V bc E^(-bc y) (θci + bc Integrate[E^(bc s) θw[x, s], {s, 0, y} ]) == 0

bcx = {D[θw[x, y], x] == 0 /. x -> 0, D[θw[x, y], x] == 0 /. x -> 1}

bcy = {D[θw[x, y], y] == 0 /. y -> 0, D[θw[x, y], y] == 0 /. y -> 1}

Can eq. $4$ subjected to the given boundary conditions, be solved using the method of finite Fourier or any other integral transform ?

Some parameter values are:

bh=0.433;bc=0.433;λh = 2.33 10^-6; λc = 2.33 10^-6; V = 1;θhi = 1; θci = 0;

NOTE: In case the above b.c. causes problems, a solution with $frac{partial theta_w(1,y)}{partial x}=frac{partial theta_w(x,1)}{partial y}=0, theta_w(x,0)=theta_{hi}, theta_w(0,y)=theta_{ci}$ is also acceptable.

---

NOTE 2: I have observed that $(4)$ can be variable separated if the ansatz $theta_w=e^{-b_h x}f(x)e^{-b_c y}g(y)$ is used. But this path lead me to nowhere. On applying SOV to the eqn. above 4 one gets:
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}$ where $F(x) := int f(x) , mathrm{d}x$ and $G(y) := int g(y) , mathrm{d}y$
with bc transformed as
with boundary conditions as:

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

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

The non-homogeneous condition in $F$ is: $beta_h e^{-beta_c y}G'(y)F(0)=1$ . I could not proceed anymore.