I'm not finding solving this enormously difficult but it is both tedious and very lengthy. For that reason I'll take the unusual approach to publish the solution in 3 parts:
Part 1: up to solving the homogeneous PDE.
Part 2: solving the non-homogeneous PDE.
Part 3: testing the solution with a simple example.
Problem statement:
$$u_t=ku_{xx}$$
With boundary and in initial conditions:
$$u(x,0)=f(x)$$
$$u(0,t)=b_1(t),u(L,t)=b_2(t)$$
So we're looking to solve the heat equation in one dimension, without heat sinks or sources but with time-dependent boundary conditions.
Part 1:
We assume the solution to be of the partitioned form:
$$u(x,t)=s(x,t)+v(x,t)$$
So that:
$$s_t+v_t=ks_{xx}+kv_{xx}...(i)$$
Translating the boundary/initial conditions:
$$s(x,0)+v(x,0)=f(x)$$
$$s(0,t)+v(0,t)=b_1(t)$$
$$s(L,t)+v(L,t)=b_2(t)$$
We assume $s(x,t)$ to be linear in $x$ but time-dependent, so the general form of $s(x,t)$ is:
$$s(x,t)=a(t)x+b(t)$$
So:
$$s_{xx}=0$$
So with $s(0,t)$ and $s(L,t)$ we get:
$$b(t)=b_1(t)$$
$$a(t)L+b(t)=b_2(t)$$
$$a(t)=\frac{b_2(t)-b_1(t)}{L}$$
So:
$$s(x,t)=\frac{b_2(t)-b_1(t)}{L}x+b_1(t)$$
With $(i)$ and $s_{xx}=0$, we get:
$$v_t=kv_{xx}-s_t$$
We solve the homogeneous equation first:
$$v_t=kv_{xx}$$
With boundary conditions:
$$v(0,t)=0,v(L,t)=0...(ii)$$
Using the Ansatz:
$$v(x,t)=X(x)T(t)$$
Obtaining separation of variables:
$$\frac1k \frac{T'}{T}=\frac{X''}{X}=-m^2$$
ODE for $X(x)$:
$$X''+m^2X=0$$
$$X=c_1\cos mx+c_2\sin mx$$
With boundary conditions $(ii)$:
$$c_1=0$$
$$m=\frac{n\pi}{L}, n=1,2,3....$$
$$X_n(x)=\sin\Big(\frac{n\pi x}{L}\Big)$$
Part 2:
Now for the non-homogeneous PDE:
$$v_t=kv_{xx}-s_t$$
With boundary/initial conditions:
$$v(0,t)=0,v(L,t)=0$$
$$v(x,0)=f(x)-s(x,0)$$
Assume the solution to be of the form:
$$v(x,t)=\displaystyle \sum_{n=1}^{+\infty}T_n(t)X_n(x)$$
$$s_t=\frac{b'_2(t)-b'_1(t)}{L}x+b'_1(t)=-\displaystyle \sum_{n=1}^{+\infty}Q_n(t)X_n(x)$$
$$Q_n(t)=-\frac{\int_0^Ls_tX_n(x)dx}{\int_0^LX_n(x)X_n(x)dx}$$
Our solution must of course obey the original equation, so:
$$\frac{\partial}{\partial t}\displaystyle \sum_{n=1}^{+\infty}T_n(t)X_n(x)=k\frac{\partial^2}{\partial t^2}\Bigg[\displaystyle \sum_{n=1}^{+\infty}T_n(t)X_n(x)\Bigg]+\displaystyle \sum_{n=1}^{+\infty}Q_n(t)X_n(x)$$
We know that:
$$X''_n(x)=-m^2X_n(x)$$
So:
$$\displaystyle \sum_{n=1}^{+\infty}T'_n(t)X_n(x)=\displaystyle \sum_{n=1}^{+\infty}-km^2T_n(t)X_n(x)+\displaystyle \sum_{n=1}^{+\infty}Q_n(t)X_n(x)$$
$$\displaystyle \sum_{n=1}^{+\infty}[T'_n(t)+km^2T_n(t)]X_n(x)]=\displaystyle \sum_{n=1}^{+\infty}Q_n(t)X_n(x)$$
So that:
$$T'_n(t)+km^2T_n(t)=Q_n(t)$$
Solve with an integration factor to:
$$T_n(t)=e^{-km^2t}\int_0^te^{km^2t}Q_n(t)dt+C_ne^{-km^2t}$$
Using the initial condition:
$$v(x,0)=f(x)-s(x,0)=\displaystyle \sum_{n=1}^{+\infty}T_n(0)X_n(x)=\displaystyle \sum_{n=1}^{+\infty}C_nX_n(x)$$
So:
$$C_n=\frac{\int_0^L[f(x)-s(x,0)]X_n(x)dx}{\int_0^LX_n(x)X_n(x)dx}$$
Part 3:
In Part 3 I'll test drive the solution for the simple case:
$$u(x,0)=u_0$$
$$u(0,t)=u_0-at$$
$$u(L,t)=u_0+bt$$
Bits and pieces we need:
$$s(x,t)=\frac{b_2(t)-b_1(t)}{L}x+b_1(t)=u_0+\frac{a+b}{L}xt-at$$
$$s_t=\frac{b'_2(t)-b'_1(t)}{L}x+b'_1(t)=\frac{a+b}{L}x-a$$
$$s(x,0)=u_0$$
$$Q_n(t)=-\frac{\int_0^Ls_tX_n(x)dx}{\int_0^LX_n(x)X_n(x)dx}$$
Because the boundary conditions are only linear in $t$ the $Q_n$ turn out time independent and compute to:
$$Q_n=2\Big[\frac{b}{n\pi}\cos n\pi+\frac{a}{n\pi}\Big]$$
Because:
$$s(x,0)=u_0=f(x)\implies C_n=0$$
So slightly reworked we get:
$$T_n(t)=\frac{2L^2}{kn^3\pi^3}\Big[(-1)^nb+a\Big]\Big[1-e^{-k\Big(\frac{n\pi}{L}\Big)^2t}\Big]$$
Put it altogether with:
$$u(x,t)=s(x,t)+\displaystyle \sum_{n=1}^{+\infty}T_n(t)X_n(x)$$
Where $X_n(x)=\sin\Big(\frac{n\pi x}{L}\Big)$
Let's evaluate for $u_0=100$, $L=1$, $k=1$, $a=1$ and $b=2$.
Firstly let's look at $s(x,t)$. The significance of $s(x,t)$ is as follows. If we allow temperature evolution of the boundaries for some time $t$, then stop heating/cooling and maintain the boundary temperatures until thermal equilibrium is achieved, then $s(x,t)$ describes the temperature profile, e.g.:
Here's $u(x,t)$ with the first three terms computed:
And here's $v(x,t)$ with the first three terms computed:
