The governing differential equation is
$$\nabla^2 T=0 \tag A$$
The boundary conditions for this problem are as foll0ws:
$$T(0,y,z)=T_{hi} \tag {1A}$$
$$T(L,y,z) = T(x,0,z) = T(x,l,z) = T(x,y,0)= T(x,y,\mu) = 0 \tag {1B}$$
I need to solve for the distribution $T(x,y,z)$. I include my attempt here.
Solution attempt
Homogeneous Dirichlet type B.C. on $y$ and $z$ faces allow us to write the following form of preliminary temperature distribution:
$$T(x,y,z) = \sum_{n,m=1}^{\infty}\sin\bigg(\frac{n\pi y}{l}\bigg)\sin\bigg(\frac{m\pi z}{\mu}\bigg)T_{nm}(x) \tag {1C}$$
We substitute $(1\mathrm{C})$ in $\mathrm{A}$ and apply the orthogonality properties of $\sin\bigg(\frac{n\pi y}{l}\bigg)$ and $\sin\bigg(\frac{m\pi z}{\mu}\bigg)$ to obtain the following:
$$\frac{\mathrm{d^2}T_{nm}(x)}{\mathrm{d}x^2} - \gamma_{nm}^2 T_{nm}(x) = 0 \tag {1D}$$ where $\gamma_{nm} = \sqrt{(\frac{n\pi}{l})^2+(\frac{m\pi}{\mu})^2}$
The general solution of $(1\mathrm{D})$ is of the form:
$$T_{nm}(x) = A_{nm}e^{\gamma x} + B_{nm}e^{-\gamma x} \tag {1E}$$
$A_{nm},B_{nm}$ are the unknown Fourier coefficients that need to be determined. Applying $T(L,y,z) = 0$ on $(1\mathrm{D})$ and using $(1\mathrm{E})$, we arrive at:
$$A_{nm}e^{\gamma L} + B_{nm}e^{-\gamma L} = 0$$
$$\bf{B_{nm} = - A_{nm}e^{2\gamma L} \tag {1F}}$$
$$T_{nm}(x) = A_{nm}(e^{\gamma x} - e^{2\gamma L - \gamma x}) \tag {1G}$$
Using the B.C. $(1\mathrm{A})$:
$$T_{hi} = \sum_{n,m=1}^{\infty}\sin\bigg(\frac{n\pi y}{l}\bigg)\sin\bigg(\frac{m\pi z}{\mu}\bigg)A_{nm}(1 - e^{2\gamma L}) \tag {1H}$$
Multiplying both sides of $(1\mathrm{H})$ by $\int_0^l\sin\bigg(\frac{k\pi y}{l}\bigg)\mathrm{d}y$ and $\int_0^{\mu}\sin\bigg(\frac{j\pi z}{\mu}\bigg)\mathrm{d}z$ and using the principle of orthogonality, we arrive at:
$$T_{hi}\int_0^l \int_0^\mu \sin\bigg(\frac{k\pi y}{l}\bigg)\sin\bigg(\frac{j\pi z}{\mu}\bigg) \mathrm{d}y\mathrm{d}z = \frac{l\mu}{4} A_{kj} (1 - e^{2\gamma L}) \tag {1I}$$
Solving the definite integrals involved , we arrive at (for any arbitrary integer $n$ and $m$):
$$\bf{A_{nm} = \frac{4 T_{hi}}{nm\pi^2(1-e^{2\gamma L}) }(1-\cos(n\pi))(1-\cos(m\pi)) \tag {1J}}$$
$$T(x,y,z) = \sum_{n,m=1}^{\infty}(A_{nm}e^{\gamma x} + B_{nm}e^{-\gamma x})\sin\bigg(\frac{n\pi y}{l}\bigg)\sin\bigg(\frac{m\pi z}{\mu}\bigg) \tag {1K}$$
On coding this solution in MATLAB. and substituting $x=0$, I find that the answer is not equal to $T_{hi}$. For parameter values $L=0.5,l=0.5,\mu=0.05,T_{hi}=50$. When i evaluate for $x=0,y=l/2,z=\mu /2$, the answer should be $T=50$, but it evaluates to $96$ using the first four terms in the series. The series surely converges.
Is there something wrong, in the way, I am doing the problem ?Any help is greatly appreciated.
An observation
When I take $z=\mu /2, y=l/2$ and $x=0$, and consider only odd values for $n$ and $m$, the solution can be written as:
$$T\bigg(0,\frac{l}{2},\frac{\mu}{2}\bigg) = \frac{16 T_{hi}}{\pi^2}\underbrace{\bigg[1 + \frac{1}{9} + \frac{1}{25} + \frac{1}{49} + ........\bigg]}_{\pi^2 /8} = 2T_{hi} \tag {1L}$$
So, is there something wrong with my analytical solution ?
Are you sure you typed it in correctly? I checked your solution against mine and they're the same. I'll leave my answer below in case it helps anyways.
Often times I find it better to start from the beginning when I'm stuck to remove doubt of errors made in assumptions. \begin{cases} \begin{align} \Delta u(x,y,z) &= 0 \\ u(0,y,z) &= u_0 \\ u\vert_{\partial\Omega\backslash\{x=0\}} &= 0. \end{align} \end{cases} You know about separation of variables already, so we may take the eigenfunctions in the $y$ and $z$ directions as $\sin \frac{n\pi y}{\ell}$ and $\sin \frac{m\pi z}{\mu}$. The equation in $x$ is known as $a_{mn}e^{\gamma_{mn}x} + b_{mn}e^{-\gamma_{mn}x}.$ Thus $u$ is the resulting infinite linear combination $$ u(x,y,z) = \sum_{m,n \geq 1}\left(a_{mn}e^{\gamma_{mn}x} + b_{mn}e^{-\gamma_{mn}x}\right)\sin\frac{n\pi y}{\ell}\sin\frac{m\pi z}{\mu}. $$ The boundary conditions on $y$ and $z$ are already satisfied so we look towards those for $x$. \begin{cases} \begin{align} u_0 &= \sum_{m,n \geq 1}(a_{mn} + b_{mn})\sin\frac{n\pi y}{\ell}\sin\frac{m\pi z}{\mu} \\ 0 &= \sum_{m,n \geq 1} \left(a_{mn}e^{\gamma_{mn}L} + b_{mn}e^{-\gamma_{mn}L}\right)\sin\frac{n\pi y}{\ell}\sin\frac{m\pi z}{\mu} \end{align} \end{cases} Note that using orthogonality is not multiplying each side by the respective integrals; it is multiplying each side by the eigenfunctions with dummy indices and then integrating. Continuing with this in mind, \begin{cases} \begin{align} \int_0^\ell \int_0^\mu u_0 \sin\frac{n\pi y}{\ell} \sin\frac{m\pi z}{\mu} \, \mathrm{d}z \,\mathrm{d}y &= \frac{\mu\ell}{4}(a_{mn} + b_{mn}) \\ 0 &= a_{mn}e^{\gamma_{mn}L} + b_{mn}e^{-\gamma_{mn}L} \end{align} \end{cases} The first integral is easily found and we have the system of equations $$ \begin{pmatrix} 1 & 1 \\ e^{\gamma_{mn}L} & e^{-\gamma_{mn}L} \end{pmatrix} \begin{pmatrix} a_{mn} \\ b_{mn} \end{pmatrix} = \begin{pmatrix} \frac{4\big((-1)^m - 1\big) \big((-1)^n - 1\big)}{\pi^2mn}u_0 \\ 0 \end{pmatrix}. $$ Inverting to solve for $a_{mn}$ and $b_{mn}$ then, \begin{align} \begin{pmatrix} a_{mn} \\ b_{mn} \end{pmatrix} &= -\frac{\operatorname{csch}(\gamma_{mn}L)}{2} \begin{pmatrix} e^{-\gamma_{mn}L} & -1 \\ -e^{\gamma_{mn}L} & 1 \end{pmatrix} \begin{pmatrix} \frac{4\big((-1)^m - 1\big) \big((-1)^n - 1\big)}{\pi^2mn} \\ 0 \end{pmatrix} \\ &= -\frac{2u_0}{\pi^2}\frac{\big((-1)^m - 1\big) \big((-1)^n - 1\big)}{\sinh(\gamma_{mn}L) mn} \begin{pmatrix} e^{-\gamma_{mn}L} \\ -e^{\gamma_{mn}L} \end{pmatrix}. \end{align} Thus $u$ is written \begin{equation} u(x,y,z) = -\frac{4u_0}{\pi^2} \sum_{m,n \geq 1} \Gamma_{mn} \frac{\sinh\big(\gamma_{mn}(x-L)\big)}{\sinh(\gamma_{mn}L)} \sin\frac{n \pi y}{\ell} \sin \frac{m\pi z}{\mu} \end{equation} with $$ \Gamma_{mn} = \frac{\big((-1)^m - 1\big) \big((-1)^n - 1\big)}{mn}. $$
$\hskip 1 in$
Mathematica code - it should be copy-pastable - enjoy.