Time series model derived from differential equation

I am trying fit a Bayesian time series model that is derived by discretizing an ordinary differential equation. Here’s a description of a simplified version of the model.

Consider an insulated box at a (time-dependent) temperature T_c in an ambient environment at a (time-dependent) temperature T_h. Assume a simple heat conduction model where the thermal resistance of the box is R and the heat capacity of the air in the box is C. Then the differential equation describing the box temperature is the following:

\frac{dT_c}{dt} = \frac{1}{RC} \left( T_h - T_c \right)

After doing a first-order discretization of the derivative, we arrive at the following difference equation:

T_c^{(k)} = T_c^{(k-1)} + \frac{\Delta t}{RC} \left( T_h^{(k-1)} - T_c^{(k-1)} \right)

Or, after collecting the T_c^{(k-1)} terms:

T_c^{(k)} = \frac{\Delta t}{RC} T_h^{(k-1)} + \left( 1 - \frac{\Delta t}{RC} \right) T_c^{(k-1)}

Now if we turn this into a Bayesian model, it looks something like the following:

R \sim N(\ldots) \\ C \sim N(\ldots) \\ \mu^{(k)} = \frac{\Delta t}{RC} T_h^{(k-1)} + \left( 1 - \frac{\Delta t}{RC} \right) T_c^{(k-1)} \\ T_c^{(k)} \sim N(\mu^{(k)}, \sigma^2)

From perusing similar threads on Discourse, it sounds like creating the T_c^{(k)} variables in a for loop is not recommended. And it appears that pm.AR cannot be used here, because the \frac{\Delta t}{RC} T_h^{(k-1)} term is not constant—it changes at each time step since T_h is time-dependent. So, is there a way to represent this model in pymc?

For recursive models, the tool in PyMC is a scan. It should be able to handle your model just fine, but there are some small details to pay attention to. Examples are here, here, and here, with some small discussion of the finer points here.

One question about the model: what is the law of motion for T_h?