Let's back up a little bit and provide a comprehensive answer to these types of problems.
Suppose $u(x,t)$ solves
\begin{align}
u_t&=u_{xx}, \qquad 0 < x < \ell,\ t>0,\\
u(0,t)&=f(t),\\
u(\ell, t)&=g(t),\\
u(x,0)&=h(x).
\end{align}
In the subsequent work, we will impose whatever smoothness conditions on the initial and boundary data we need to get convergence of the involved series.
First, standard separation of variables shows that the solution to the problem with homogeneous BCs is $$u(x,t)=\sum_{n=1}^\infty b_n\sin(\sqrt{\lambda_n}\,x)e^{-\lambda_n t},$$ where $\lambda_n=(n\pi/\ell)^2$, $n=1,2,\dots$ In other words, for each fixed $t>0$,
$$u(x,t)=\sum_{n=1}^\infty u_n(t)\sin(\sqrt{\lambda_n}\,x)\quad\text{where}\quad u_n(t)={2\over \ell}\int_0^\ell u(x,t)\sin(\sqrt{\lambda_n}\,x)\,dx.$$ (This is a key observation. I hope the notation isn't confusing: $u_n$ represents the coefficients in the series for $u$, not a partial derivative.)
Then, differentiating the series above, define $v_n(t)$ and $w_n(t)$ as the coefficients in the series for ${\partial u\over \partial t}$ and ${\partial^2 u\over \partial x^2}$, respectively:
\begin{align}
{\partial u\over \partial t}&=\sum_{n=1}^\infty v_n(t)\sin(\sqrt{\lambda_n}\,x)
\quad\text{where}\quad
v_n={2\over \ell}\int_0^\ell {\partial u\over \partial t}\sin(\sqrt{\lambda_n}\,x)\,dx={du_n\over dt},\\
{\partial^2 u\over \partial x^2}&=\sum_{n=1}^\infty w_n(t)\sin(\sqrt{\lambda_n}\,x)
\quad\text{where}\quad
w_n={2\over \ell}\int_0^\ell {\partial^2 u\over \partial x^2}\sin(\sqrt{\lambda_n}\,x)\,dx.
\end{align}
Integrating the $w_n(t)$ by parts, simplifying the trig terms, and applying the BCs, we get
\begin{align}
w_n(t)&=-{2\over \ell}\int_0^{\ell} \lambda_n u(x,t)\sin(\sqrt{\lambda_n}\,x)\,dx\\
&\qquad\qquad+{2\over \ell}\left[u_x(x,t)\sin(\sqrt{\lambda_n}\,x)-\sqrt{\lambda_n}\,u(x,t)\cos(\sqrt{\lambda_n}\,x)\right]\Bigg|_{x=0}^{x=\ell}\\
&=-\lambda_nu_n(t)+\underbrace{{2\sqrt{\lambda_n}\over \ell}\left[f(t)+(-1)^{n+1}g(t)\right]}_{F(t)}.
\end{align}
From the PDE, $$u_t=u_{xx}\implies v_n(t)=w_n(t)\implies {du_n\over dt}=-\lambda_n u_n(t)+F(t),$$ and thus the coefficients $u_n(t)$ we seek are found by solving the (ODE!) initial-value problem
\begin{align}
{du_n\over dt}+\lambda_nu_n(t)&=F(t),\\
u_n(0)&={2\over \ell}\int_0^\ell h(x)\sin(\sqrt{\lambda_n}\,x)\,dx,
\end{align}
by the method of your choice.
This is called the method of eigenfunction expansions. Transform methods are also available, but that is a separate post.
Hope that helps.
I'm sorry, but I cannot tell where you went wrong; your answer makes no sense from a dimensional standpoint. As for the convergence of the series, I will answer below, but first let's actually solve the equation correctly.
EXACT SOLUTION OF THE EQUATION
In my humble opinion, I think these particular mixed BCs do not lend themselves easily to a simple solution by separation of variables (although the final solution is indeed of that form). I chose, rather, to solve this equation using a Laplace transform over $t$; that is
$$\hat{u}(x,s) = \int_0^{\infty} dt \, u(x,t) e^{-s t}$$
Then the PDE becomes on ODE on $x$:
$$\alpha \frac{d^2}{d x^2} \hat{u}(x,s) - s \hat{u}(x,s)=0$$
$$\left[\frac{d}{dx} \hat{u}(x,s)\right]_{x=0} = -\frac{1}{s}$$
$$\hat{u}(L,s) = 0$$
Note that the ODE already incorporates the initial condition.
The solution to this equation (the details of which solution I leave to the reader) is
$$\hat{u}(x,s) = \sqrt{\alpha} \frac{\sinh{\left[\frac{L-x}{\sqrt{\alpha}} \sqrt{s}\right]}}{s^{3/2} \cosh{\left[\frac{L}{\sqrt{\alpha}} \sqrt{s}\right]}}$$
Getting the solution to the original equation is then a matter of taking an inverse Laplace transform:
$$u(x,t) = \frac{1}{i 2 \pi} \int_{c-i \infty}^{c+i \infty} ds \, e^{s t} \hat{u}(x,s)$$
Although the integration contour appears to be complicated by the branch points in the LT expression, the ILT turns out to be simply equal to the sum of the residues of the poles of $\hat{u}(x,s)$. The reason this is so is because the functions of $\sqrt{s}$ are all even, so that there is no effect of crossing any branch cuts. The poles of $\hat{u}$ are at $s = -(n+1/2)^2 \pi^2 \alpha/L^2$ for integers $n$ and at $s=0$. Using the fact that, for a function of the form $f(s) = p(s)/q(s)$, the residue of a simple pole $s_0$ of $f$ is $p(s_0)/q'(s_0)$, I get for the ILT, after some algebra:
$$u(x,t) =L-x- 2 \frac{L}{\pi^2} \sum_{n=0}^{\infty} \frac{1}{(n+1/2)^2} \cos{\left[ \left ( n+\frac12 \right) \pi \frac{x}{L}\right]} e^{-\alpha (n+1/2)^2 \pi^2 t/L^2}$$
The series obviously converges, and you can perform some rudimentary checks to see that the initial and boundary conditions are satisfied, e.g., note that
$$\sum_{n=0}^{\infty} \frac{1}{(n+1/2)^2} = \frac{\pi^2}{2}$$
Also, you should check that the solution jibes with intuition at the steady state (i.e., as $t \to \infty$): here, $u(x,t) \sim L-x$, which makes physical sense (the temperature gradient at $x=0$ implies heat flow out of the left end of the bar).
RATE OF CONVERGENCE
The series above, while clearly convergent, does so painfully slow, especially for small $t$. Even for $n=50$, there is significant oscillation in the sum for $t=0$ which should ultimately converge to zero everywhere in $(0,L)$. Nevertheless, as $t$ increases, I do see increasingly fast convergence. The following is the behavior of the solution for $L=\alpha=1$, $n=10$ terms for times between $t=0.01$ and $t=5$:
![enter image description here](https://i.stack.imgur.com/r8XKf.jpg)
$n=20$ is identical.
Best Answer
The equation in $Z$ is $$ -Z'' = \lambda^2 Z, \\ Z'(0)-EZ(0)=0 \\ Z'(1)+DZ(1)=0. $$ If $Z_1$ is a solution for $\lambda_1$ and $Z_2$ is a solution for $\lambda_2$, then \begin{align} (\lambda_2^2-\lambda_1^2)\int_{0}^{1}Z_1(z)Z_2(z)dz & = \int_{0}^{1}Z_1Z_2''-Z_1''Z_2 dz \\ & = \int_{0}^{1}\frac{d}{dz}(Z_1Z_2'-Z_1'Z_2)dz \\ & = Z_1Z_2'-Z_1'Z_2|_{0}^{1} \\ & = \left.\left|\begin{array}{cc}Z_1 & Z_2 \\ Z_1' & Z_2'\end{array}\right|\right|_{0}^{1} = 0. \end{align} The determinants at $0$ and at $1$ are individually $0$ because the matrices have non-trivial null spaces: $$ \left[\begin{array}{cc} Z_1(0) & Z_1'(0) \\ Z_2(0) & Z_2'(0)\end{array}\right]\left[\begin{array}{c}1 \\ -E\end{array}\right] = 0, \\ \left[\begin{array}{cc} Z_1(1) & Z_1'(1) \\ Z_2(1) & Z_2'(1)\end{array}\right]\left[\begin{array}{c}1 \\ \;\;D\end{array}\right] = 0. $$ Therefore, if $\lambda_1\ne\lambda_2$, $$ \int_{0}^{1}Z_1(z)Z_2(z)dz = 0. $$