Smells of homework, so just a hint :
The problem here is that we have a "nuisance" unknown function : the temperature, which is a function of the time. in order to be able to solve (numerically), we must be able to compute this temperature for any point in $\left[0~24\right]$. This can be accomplished vie the use of a interpolating function such as spline()
This Sagecell example illustrate the result, which can be used to build the call to rk4 :
T=[21,20,17,16,18,21,25,27,28,26,23,21,21]
d=[0,2,4,6,8,10,12,14,16,18,20,22,24] # or d=[2*u for u in range(len(T))]
P=zip(d,T)
S=spline(P)
plot(lambda t:S(t),(0,24))+points(P)
EDIT : You should state C(0) in order to get a numerical solution...
EDIT 2: If you are able and willing to state that your temperature data are periodic of period 24, it might be preferable to create an interpolation function based on the FFT of your data. It should be possible (and efficient) to create such a "periodic spline" in Sage, but, being lazy, I'll illustrate how to use R's splinefun function in this Sagecell example :
T=[21,20,17,16,18,21,25,27,28,26,23,21,21]
d=[0,2,4,6,8,10,12,14,16,18,20,22,24] # or d=[2*u for u in range(len(T))]
rCall="splinefun(x=c({}), y=c({}), method='periodic')"
rCall=rCall.format(", ".join([str(u) for u in d]),
", ".join([str(u) for u in T]))
foo=r(rCall)
plot(lambda u:foo(u).sage(),(-24,48)) + points(zip(d,T)) + points(zip([u-24 for u in d],T)) + points(zip([u+24 for u in d],T))
I'm not new to Sage but since I don't use it a lot I'm kinda an inexperienced user. But what I would do (whether on Sage or any other numerical package) is to find Temp(t), that is, find Temp as a function of 't' given the 't' and 'Temp' points you supplied as lists. I plotted your points and they represent a broken sinusoidal line, but there may be ways to find a reasonable expression for Temp(t). Then you would substitute Temp(t) in your code by the expression you found.