On the last step of solving a three-dimensional Laplace equation,($\nabla^2T=0$) with BC(s) as $T(0,y,z) = T(L,y,z) = T_a$,
$T(x,0,z) = T(x,l,z) = T_a$,
$\frac{\partial T(x,y,0)}{\partial z} = p_c\bigg(e^{-b_cy/l}\left[T_{ci} + \frac{b_c}{l}\int_0^y e^{b_cs/l}T(x,s,z)ds\right] - T(x,y,0)\bigg)$
$\frac{\partial T(x,y,w)}{\partial z} = p_h\bigg(e^{-b_hx/L}\left[T_{hi} + \frac{b_h}{l}\int_0^x e^{b_hs/l}T(s,y,z)ds\right] - T(x,y,w)\bigg)$
On using a solution form
$T(x,y,z)=\sum_{m,n=1}^{\infty}(A_{nm}e^{\gamma z} + B_{nm}e^{-\gamma z})\sin(\frac{n\pi x}{L})\sin(\frac{m\pi y}{l})+T_a.$
I need to evaluate the Fourier Coefficients $A_{nm}$ and $B_{nm}$ using the following equation which was obtained after applying the $z=0$ BC
$$ \frac{1}{p_c}\sum_{n,m=1}^\infty\sin(\frac{n\pi x}{L})\sin(\frac{m\pi y}{l})\gamma(A_{nm}-B_{nm}) = \delta e^{\frac{-b_c y}{l}} + P-Q+R-S \rightarrow (\mathrm{1}) $$
$P =\sum_{n,m=1}^\infty\frac{b_c^2}{b_c^2+(m\pi)^2}(A_{nm}+B_{nm})\sin(\frac{n\pi x}{L})\sin(\frac{m\pi y}{l})$
$Q =\sum_{n,m=1}^\infty \frac{b_c m\pi}{b_c^2+(m\pi)^2}(A_{nm}+B_{nm})\sin(\frac{n\pi x}{L})\cos(\frac{m\pi y}{l})$
$R = \sum_{n,m=1}^\infty \frac{b_c m\pi}{b_c^2+(m\pi)^2}(A_{nm}+B_{nm})e^{\frac{-b_c y}{l}}$
$S = \sum_{n,m=1}^\infty(A_{nm}+B_{nm})\sin(\frac{n\pi x}{L})\sin(\frac{m\pi y}{l})$
Here ,$\gamma^2 = (\frac{n\pi}{L})^2 + (\frac{m\pi}{l})^2$ and $c,\delta=(T_{ci}-T_a)$ are constants.
Obviously, there would be another relation that would be a linear equation in $A$ and $B$ to finally deduce the values of $A_{nm}$ and $B_{nm}$ separately.My question pertains to the fact as to how should I handle the $e^{\frac{-b_c y}{l}}$ here.
Attempt
I can multiply both sides of the equation $\mathrm{(1)}$ with $\int_0^L\sin(\frac{k\pi x}{L})$ and $\int_0^l\sin(\frac{j\pi y}{l})$ and then use the principle of orthogonality to finally arrive at
For some $n=k$ and $m=j$, $\mathrm{(1)}$ becomes
$$ \frac{1}{p_c}\frac{L}{2}\frac{l}{2}\gamma(A_{kj}-B_{kj})=\delta e^{\frac{-b_c y}{l}}\int_0^L\sin(\frac{k\pi x}{L})\mathrm{d}x\int_0^l\sin(\frac{j\pi y}{l})\mathrm{d}y+\underbrace{\frac{b_c^2}{b_c^2+(j\pi)^2}(A_{kj}+B_{kj})\frac{L}{2}\frac{l}{2}}_{P} - \overbrace{0}^{Q} + \underbrace{\sum_{n,m=1}^\infty \frac{b_c m\pi}{b_c^2+(m\pi)^2}(A_{nm}+B_{nm})e^{\frac{-b_c y}{l}}\int_0^L\sin(\frac{k\pi x}{L})\mathrm{d}x\int_0^l\sin(\frac{j\pi y}{l})\mathrm{d}y}_{R} - \overbrace{(A_{kj}+B_{kj})\frac{L}{2}\frac{l}{2}}^{S} $$
Fourier expansion of $e^{\frac{-b_c y}{l}}$ in the interval $y \in [0,l]$
I expanded $e^{\frac{-b_c y}{l}}$
$$e^{\frac{-b_c y}{l}} = \frac{(1-e^{-b_c})}{b_c} + \sum_{r=1}^{\infty}\bigg[\frac{2b_c(1-e^{-b_c})}{(b_c)^2 + (2r\pi)^2}\cos\bigg(\frac{2r\pi y}{l}\bigg) + \frac{4r\pi(1-e^{-b_c})}{(b_c)^2 + (2r\pi)^2}\sin\bigg(\frac{2r\pi y}{l}\bigg)\bigg]$$
I guess, I need to substitute this in $(1)$ for $e^{\frac{-b_c y}{l}}$
Any suggestions would be helpful and really appreciated.
I would mention that I am particularly looking for some guiding points or a discussion, in case a full answer is not possible.
