It might be that you would like to look into this book
https://www.amazon.de/Conduction-Solids-Oxford-Science-Publications/dp/0198533683
You can achieve the given solution by the method of Laplacetransformation. This gives an ODE instead of the PDE for heat-transport. I will try to line it out.
Write the heat equation ( PDE )as
heateq = D[fT[x, t], t] - a D[fT[x, t], x, x]
and do a Laplace Transform
LaplaceTransform[heateq, t, s]
Now note that LaplaceTransform[fT[x, t], t, s] is a function of x only, we have an ordinary differential equation. Substituting LaplaceTransform[fT[x, t], t, s] with fTL[x] rewrite the ODE and solve for fTL
sol = DSolve[-T0 + s fTL[x] - a D[fTL[x], x, x] == 0, fTL[x], x] // Flatten
and get
fT1 = fTL[x] /. sol
fT1 has to be bounded (not grow above all limits for growing x ), so C[1] has to be 0
fT2 = fT1 /. C[1] -> 0
and for x = 0 your solution has to be T1 ( in the L-regions this is T1/s ),
cond = (fT2 /. x -> 0) == T1/s
this gives the value for C[2] and the final result (in L-space)
fT3 = fT2 /. (Solve[cond, C[2]] // Flatten)
Now transform your result back into the time-domain (which may be next to impossible)
fxt = Simplify[InverseLaplaceTransform[fT3, s, t], x > 0 && a > 0]
and get the result you asked for.