Suppose that $u$ is the solution to the heat equation with mixed Neumann and Robin boundary conditions
\begin{align}
&\partial_tu(t,x) = k \partial_{xx} u(t,x), &&t>0, 0<x<L, \\
&\partial_x u(t,0) = \phi(t), &&t>0,\\
&\partial_x u(t, L)= \alpha u(t,L), &&t>0, \\
&u(0, x)= u_0(x), &&0<x<L,
\end{align}
where $\phi$ is a smooth function and $k > 0$. If $\phi(t) = H(t-t_0)$, i.e., $\phi(t) = 1$ for $t \geq t_0$ and $0$ otherwise, and disregarding the initial condition for a second, I can find a solution to the equation which - if I'm not mistaken - reads for all $t > t_0$
$$
u(t, x; t_0) = u_s(x)+\sum_{n\geq 1} c_ne^{-k\mu_n^2(t-t_0)}\cos(\mu_nx),
$$
for appropriate coefficients $c_n$, where $\mu_n$, $n =1, 2, \ldots$, are the solutions to
$$
-\frac{\alpha}{\mu}=\tan(\mu L),
$$
and where the stationary solution $u_s$ is given by
$$
u_s(x) = x- L + \frac1\alpha.
$$
I'm tempted to then use a superposition principle and write (formally)
$$
\phi(t) = \int \delta(t - t') \phi(t') dt' = \int H(t-t') \phi'(t') dt',
$$
so that
\begin{align}
u(t,x) &= \int u(t, x;t') \phi'(t') dt' \\
&= u_s(x) \int \phi'(t')dt' + \sum_{n\geq 1} c_n \cos(\mu_nx) \int e^{-k\mu_n^2(t-t')} \phi'(t') dt',
\end{align}
The coefficients $c_n$ can then be found by imposing the initial condition to be satisfied.
I left the extrema of integration of the integrals above on purpose blank, because I'm a bit rusty on distribution theory and would need some help in making the above rigorous. How should I proceed in this situation?