Three dimensional Laplace equation with constant Temp. on one face. [Solution not satisfying BC] 
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 ?
 A: 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.
ClearAll["Global`*"]
\[Gamma]mn = Sqrt[((n \[Pi])/\[ScriptL])^2 + ((m \[Pi])/\[Mu])^2];
(* u0 not specified yet *)
L = 1;
\[ScriptL] = 1;
\[Mu] = 1;
(* Keep q low (~10-15 or so) as the number of terms grows as q^2 *)
(* i.e. higher q gives more accuracy but takes longer to run in effect *)
q = 5;

(*The coefficients amn and bmn {amn,bmn}={a,b}/.Solve[{a+b\[Equal]4/\
\[Pi]^2(((-1)^m-1)((-1)^n-1))/(m n)u0,a Exp[\[Gamma]mn L]+b Exp[-\
\[Gamma]mn L]\[Equal]0},{a,b}][[1]];*)

u[x_, y_, 
   z_] = (-4 u0)/\[Pi]^2 Sum[
    Sum[(((-1)^m - 1) ((-1)^n - 1))/(m n) Sinh[\[Gamma]mn (x - L)]/
       Sinh[\[Gamma]mn L] Sin[n \[Pi] y/\[ScriptL]] Sin[
       m \[Pi] z/\[Mu]], {n, 1, q}], {m, 1, q}];

Plot3D[u[0, y, z]/u0, {y, 0, \[ScriptL]}, {z, 0, \[Mu]}, 
 PlotLabel -> 
  "Laplace's Equation Solution\nat x=0 with\n\[ScriptL]=1, \[Mu]=1, \
25 terms", 
 AxesLabel -> {"y", "z", "\!\(\*FractionBox[\(u\), \(u0\)]\)"}, 
 Boxed -> False]

(* Something extra - it's computationally expensive (3d), at least for my \
computer,  so let it run for a bit *)
Manipulate[
 ContourPlot[u[x, y, z]/u0, {x, 0, L}, {y, 0, \[ScriptL]}, 
  ColorFunction -> "DarkRainbow"], {z, 0, \[Mu]}]

(* alt+.to quit if it takes too long & reduce q *)
ContourPlot3D[
 u[x, y, z]/u0, {x, 0, L}, {y, 0, \[ScriptL]}, {z, 0, \[Mu]}, 
 Mesh -> False, ColorFunction -> "Rainbow", 
 PlotLabel -> "Contours of Temperature", AxesLabel -> {"x", "y", "z"}]

