Evaluating Fourier coefficients when Exponential terms are present

I am trying an analytical solution for a 3D Laplace equaition, and just wanted to know if I could use the symbolic calculation capabilities of Mathematica to calculate the Fourier coefficients

Context (Please bear with the length)

On the last step of solving a three-dimensional Laplace equation,($$\nabla^2T=0$$, where $$\nabla^2 = \frac{\partial^2}{\partial x^2} + \frac{\partial^2}{\partial y^2} + \frac{\partial^2}{\partial z^2}$$) 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}$$ (from the $$z=w$$ BC) 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}}$$

If you look at $$(1)$$, you can see how if the expansion of $$e^{\frac{-b_c y}{l}}$$ is substituted the $$R$$ term would go to a triple summation which makes it way complicated. Can there be some alternative way to handle these calculations in Mathematica. I do know that the PDE solvers could give a result using discretization schemes. I was hoping if I could use Mathematica to help me in the final step of finding these Fourier coefficients.