Heat PDE with unusual, inhomogeneous Neumann boundary condition

The steady state solution for the original problem is $y_E(x)=\frac{\alpha \beta}{1-\alpha L} x$. The transient solution is given by $z(x,t)=y(x,t)-y_E(x)$ which now solves the PDE, for $(x,t)\in (0, L)\times (0,\infty)$, $$z_t=k z_{xx},$$ and with BCs $$z(0,t)=0 \,\text{ and } \, z_x(L,t)-\alpha z(L,t)=0$$ for $t>0$ and IC $z(x,0)=g(x)-y_E(x)$ (where $g$ is the IC of the original problem, unspecified in the OP). Thus the $\beta$ term has vanished as indicated in my comment on the OP. The transient solution can then be found by separation of variables.

For then we get two ODEs one in $x$, $\phi''+\lambda^2 \phi = 0$ for $0<x<L$ with BC $\phi(0)=0$ and $\phi'(L)-\alpha \phi(L)=0$ and one in $t$, $T'+\lambda^2 k T=0$ with the IC. Solving the first ODE and imposing the first BC gives $\phi = c_2 \sin(\lambda x),$ and imposing the second BC and avoiding trivial solutions requires $\lambda$ to solve $\tan(\lambda L)=\lambda/\alpha$ which has infinite solutions $\lambda_n$ for $n\geq 1$. All together, we get $$z(x,t)=\sum_{n} b_n \sin(\lambda_n x)e^{-\lambda_n^2 k t},$$ with initial condition $$z(x,0)=\sum_n b_n \sin(\lambda_n x)=g(x)-y_E(x),$$ which leads to $$b_n=\frac{\int_0^L [g(x)-y_E(x)]\sin(\lambda_n x) dx}{\int_0^L \sin^2(\lambda_n x) dx},$$ so that finally, returning back to $y=z+y_E$, we have $$y(x,t)=\frac{\alpha \beta}{1-\alpha L} x + \sum_n b_n \sin(\lambda_n x) e^{-\lambda_n^2 k t},$$ where $\lambda_n$ and $b_n$ are defined above.

---

Please comment for clarifications or corrections.