Nearly a year ago, I gave in the comments a non-answer that would justify calling the quantity
$$
\int u(x,t)^2~dx = E(t)
$$
the energy of the rod. Namely, one can think of minimizing the energy functional as satisfying a conservation law, and thus minimizing the energy leads you to solutions of the associated partial differential equations.
One year out, I think I have a more concrete, dimensional-analysis based answer. Let us examine the equation
$$
\partial_t u = -\alpha\partial_{xx}u
$$
for its dimensions. $u = u(x,t)$ represents temperature at the spacetime coordinate $(x,t)$, so it has units of temperature $T$. $\alpha$ here is the thermal diffusivity, which has units of length squared over time, $L^2/\tau$. Finally, not in the equation but lurking in the background is the kinetic energy, which is related to temperature by the equation
$$
E_k = \frac{3}{2}kT,
$$
where $k$ is the Boltzmann constant and has units of energy over temperature, $E/T$.
Let me use $[unit]$ to represent the dimensions of any units in any expression. So for example, the units in the heat equation check out: $\partial_t u$ has units of temperature over time, while $\partial_xx u$ has units of temperature over length squared, which gives us
$$
[T\tau^{-1}] = [L^2\tau^{-1}][TL^{-2}].
$$
Now, in our scenario the total length of the rod is fixed, so without loss of generality we may treat it as a constant. Dimensionally, this is saying we will take $[L] = [1]$, so it drops out of any of our equations. This is all fine, as we see in the heat equation the units of length cancel each other anyway.
So let us turn to the energy functional. In terms of units,
$$
\int u(x,t)^2~dx = \int [T^2]~d[L] = [T^2L] = [T^2].
$$
So the energy functional has units $[T^2]$, which isn't where we want to be (we want units of energy). However, the primary thing we do with taking the energy functional is minimize it. Minimizing the energy functional $E(t)$ is equivalent to minimizing the square root of the energy functional, $\sqrt{E(t)}$, which has units of $[T]$.
Still not there. But aha! Energy and temperature, $[E]$ and $[T]$, are related by a proportionality law! That is,the equation
$$
E_k = \frac{3}{2}kT
$$
tells us precisely how to convert from temperature to energy, and the conversion preserves order. So if we can minimize the square root of the "energy functional" $\sqrt{E(t)}$, which has units of $T$, then we automatically know how to minimize the actual energy, which according to our dimensional analysis must be something like
$$
\sqrt{\int \frac{9}{4}k^2 u(x,t)^2~dx } = \frac{3k}{2}\sqrt{E(t)},
$$
where the right-hand side is actually integrating the physical notion of energy squared. And now, the units do indeed check out: these are quantities with units of energy. There are surely other ways to define a natural notion of a physical energy functional using an integral, but the units of this definition work and it has nice mathematical properties (pretty much exactly those of $\sqrt{E(t)}$, which we know has an excellent mathematical theory). This revolves around the observation that since $k$ is a constant, minimizing the "energy" functional defined using units of heat is equivalent to minimizing what you would expect to be the "physical energy" functional.
Morally, the fact that $k$ is a constant (despite having units of $[ET^{-1}]$) gives us a temperature-energy equivalence law, which tells us that for dimensional analysis purposes $[E]$ and $[T]$ are indistinguishable. This is what the physicist in me needs to say to satisfy himself. The mathematician in me chooses units for temperature and energy so that $k = 2/3$, and then I just have the simple equivalence $E_k = T$, and now I happily exchange temperature and energy at will. And presumably this is what mathematicians of the past have done.
Treating $L$ like a constant also needs some justification, but I think if we think of $L$ as very small and argue that computing the energy of large homogeneous bodies consists of computing the energy on small pieces and summing, then we can justify that assumption as well.
You're on the right track. What you should have now is $\Theta_{xx} - s\Theta = -T_0$ (note: you have a sign error). For each fixed $s$ this is a constant-coefficient second-order linear ODE in $x$. The general solution is given by the sum of the general solution to the homogeneous equation and any particular solution to the inhomogeneous equation. There shouldn't be a need to consider $s < 0$, as the Laplace variable is usually assumed $>0$ by definition. So for $s>0$ the solution to the homogeneous equation is given as a sum of exponentials $c_1e^{\sqrt{s}x} + c_2e^{-\sqrt{s}x}$ and by inspection $\Theta_p(x) = \frac{T_0}{s}$ is a solution to the inhomogeneous equation. You can use the initial-value theorem for the Laplace transform ($f(0^+) = \lim_{s\to\infty} sF(s)$) to show that $c_1 = 0$. The boundary condition $\theta(0,t) = 0$ implies $\Theta(0,s) = 0$ for all $s>0$, which then implies $c_2 = -\frac{T_0}{s}$. Altogether from here we obtain the Laplace transform
$$ \Theta(x,s) = -T_0\frac{e^{-\sqrt{s}x}}{s} + \frac{T_0}{s}.$$
We then invert this Laplace transform. The second term just gives a unit step function, and while the inverse Laplace transform of the first term can't be expressed in terms of elementary functions, we can express it using the rule $$F(s)/s = \mathcal{L}\left[\int_0^t f(\tau)~d\tau\right](s)$$
which gives
$$
\theta(x,t) = -T_0\int_0^t \mathcal{L}^{-1}[e^{-\sqrt{s}x}](\tau)~d\tau + T_0u(t),
$$
where $u(t)$ is the unit step function. The inverse Laplace transform in the integral is a bit messy: Wolfram gives
$$
\mathcal{L}^{-1}[e^{-\sqrt{s}x}](\tau) = \frac{xe^{-\frac{x^2}{4\tau}}}{2\sqrt{\pi}\tau^{3/2}}.
$$
Making the change of variables $\eta = x/2\sqrt{\tau}$, or equivalently $\tau = x^2/4\eta^2$, $d\tau = -\frac{x^2}{2\eta^3}d\eta$ gives
$$
\int_0^t \mathcal{L}^{-1}[e^{-\sqrt{s}x}](\tau)~d\tau = \frac{2}{\sqrt{\pi}}\int_{x/2\sqrt{t}}^\infty e^{-\eta^2} ~d\eta = \operatorname{erfc}(x/2\sqrt{t}),
$$
where $\operatorname{erfc}$ denotes the complementary error function. Therefore the solution comes out to
$$
\theta(x,t) = -T_0\operatorname{erfc}(x/2\sqrt{t}) + T_0u(t).
$$
To check that this is consistent with the initial and boundary conditions: at $t=0$, $\operatorname{erfc}(x/\sqrt{t}) = 0$ for all $x>0$, so $\theta(x,0) = T_0$ for all $x>0$, while at $x=0$ one can explicitly calculate the value of $\operatorname{erfc}(0)$ (it is half of the famous Gaussian integral) to find that $\theta(0,t) = -T_0 + T_0 = 0$.
Best Answer
It is a consequence of the physics. The heat flow $Q$ is proportional to the gradient of the temperature $u$:
$$Q = -k \mathbf{\nabla}u$$
where $k$ is a heat conductivity. In a perfectly insulated material, there is no heat flow. Hence, the temperature gradient is zero. In 1D, this translates to the spatial derivative of $u$, as stated in the problem.