I have the 3D heat equation (Laplace equation)
$$\nabla^{(3)}T_s=0$$
where $\nabla^{(3)}=(\frac{\partial^{2}}{\partial x^2}+\frac{\partial^{2}}{\partial y^2}+\frac{\partial^{2}}{\partial z^2})$ defined on $x \in [0,1]$, $y \in [0,1]$ and $z \in [-1,0]$.
The boundary conditions are as follows
$$\frac{\partial T_s(0,y,z)}{\partial x}=\frac{\partial T_s(1,y,z)}{\partial x}=0$$
$$\frac{\partial T_s(x,0,z)}{\partial y}=\frac{\partial T_s(x,1,z)}{\partial y}=0$$
$$\frac{\partial T_s(x,y,0)}{\partial z}=\frac{b_cC_cw}{klL}\bigg(b_ce^{-b_cy}\int(e^{b_cy}T_s)\mathrm{d}y - T_s\bigg)$$
$$\frac{\partial T_s(x,y,-1)}{\partial z}=\frac{b_hC_hw}{klL}\bigg(T_s - b_he^{-b_hx}\int(e^{b_hx}T_s)\mathrm{d}x\bigg)$$
$b_h.b_c,C_h,C_c,w,k,l,L$ are all constants $>0$.
Any help or guiding points in solving this problem is appreciated. It is the integrals in the boundary conditions that is completely new to me.
Attempt
We can write $T_s$ as $T_s=X(x)Y(y)Z(z)$ to get the following
$$\frac{1}{X}\frac{d^2 X}{dx^2}=-\lambda_x ^{2}$$ $$\frac{1}{Y}\frac{d^2 Y}{dy^2}= \lambda_y ^{2}$$ $$\frac{1}{Z}\frac{d^2 Z}{dy^2}=-\lambda_z ^{2}$$
so that $-\lambda_x^2+\lambda_y^2-\lambda_z^2=0$
One possible combination of the Eigen functions could be hence:
$X \sim \cos(\lambda_x x)$ or $X \sim \sin(\lambda_x x)$
$Y \sim \cosh(\lambda_y y)$ or $Y \sim \sinh(\lambda_y y)$
$Z \sim \cos(\lambda_z z)$ or $Z \sim \sin(\lambda_z z)$
I cannot proceed from here on applying the boundary conditions provided for this problem.
Attempt 2
After suggestion from @Dylan and assuming a solution of the form
and assuming $p_1 = \frac{b_h C_h w}{klL}$ and $p_2 = \frac{b_c C_c w}{klL}$
$$T_w = \sum_{n,m} \cos{(n\pi x)}\cos{(m\pi y)}[A_{n,m}\cosh(\sqrt{n^2 + m^2}\pi z) + B_{n,m}\sinh(\sqrt{n^2 + m^2}\pi z) ]$$
On substituting this form in the last two boundary conditions for $z$ viz. at $z=0$ and $z=-1$ and remembering that
$\sin(n\pi)=0$, $\sinh(0)=0$, $\cosh(0)=1$, $\sinh(-x)=-\sinh(x)$ and $\cosh(-x)=\cosh(x)$, I arrive at the following two equations for $A_{m,n}$ and $B_{m,n}$
$$ B_{n,m}\bigg[p_1\bigg(\cos(n\pi x)-b_h e^{-b_h x}\bigg(\frac{b_h e^{b_h}\cos(n\pi)-b_h}{b^2 + n^2 \pi^2}\bigg)\bigg) - \pi\sqrt{n^2 + m^2}\cos(n\pi x)\bigg] + A_{n,m}p_1\bigg[\cos(n\pi x) - b_h e^{-b_h x}\bigg(\frac{b_h e^{b_h}\cos(n\pi)-b_h}{b^2 + n^2 \pi^2}\bigg)\bigg] = 0 $$
and
$$ B_{n,m}\bigg[p_2\bigg(-\sinh(\pi\sqrt{n^2 + m^2}) + b_c e^{-b_c y}\bigg(\frac{b_c e^{b_c}\cos(m\pi)-b_c}{b_c^2 + m^2 \pi^2}\bigg)-\cos(m\pi y)\bigg) - \pi\sqrt{n^2 + m^2}\cos(m\pi y)\cosh(\pi\sqrt{n^2 + m^2})\bigg] + A_{n,m} \bigg[p_2\bigg(\cosh(\pi\sqrt{n^2 + m^2}) + b_c e^{-b_c y}\bigg(\frac{b_c e^{b_c}\cos(m\pi)-b_c}{b_c^2 + m^2 \pi^2}\bigg)-\cos(m\pi y)\bigg) + \pi\sqrt{n^2 + m^2}\cos(m\pi y)\sinh(\pi\sqrt{n^2 + m^2})\bigg] = 0 $$
As can be seen we have two equations where the unknowns are $A_{n,m}$ and $B_{n,m}$. Also this is a system of homogeneous equation.
So these two must be solved simultaneously to determine $A$ and $B$ ? But won't they just give out trivial solutions ?
Note I have used $$\int_{0}^{1}e^{bx}cos(n\pi x) \mathrm{d}x = \frac{\pi e^b n \sin(\pi n) + e^{b}b\cos(\pi n) - b }{b^2 + \pi^2 m^2}$$
Attempt3
From the physical problem that this equation defines, I have concluded that there can be two other boundary conditions in place of the last two (the first four Neumann conditions remain the same)
At $x=0,y,z=-1$
$$b_h e^{-b_h x} \int e^{b_h x}T_s \mathrm{d}x = \theta_{h,i}$$
and At $x,y=0,z=0$
$$b_c e^{-b_c y} \int e^{b_c y}T_s \mathrm{d}y = \theta_{c,i}$$
Contextual information
The variable $T_s$ appears in two other equations:
$$\frac{\partial \theta_h}{\partial x} + b_h (\theta_h - T_s) = 0$$ $$\frac{\partial \theta_c}{\partial y} + b_c (\theta_c - T_s) = 0$$ with bc as
$\theta_h(0,y) = \theta_{h,i}$
$\theta_c(x,0) = \theta_{c,i}$
From here :
$\theta_h = b_h e^{-b_h x} \int e^{b_h x} T_s \, \mathrm{d}x$
$\theta_c = b_c e^{-b_c y} \int e^{b_c y} T_s \, \mathrm{d}y$
This hence explains the bc mentioned in the Attempt3
The boundary conditions for my original attempt are actually:
$$\frac{\partial T_s(x,y,0)}{\partial z}=\frac{b_cC_cw}{klL}\bigg(\theta_c - T_s\bigg)$$
$$\frac{\partial T_s(x,y,-1)}{\partial z}=\frac{b_hC_hw}{klL}\bigg(T_s - \theta_h\bigg)$$
