From a practical point of view, I can easily integrate the equations, determine the state trajectory, then call a function that evaluates

I am generating some ODE's that I integrate with scipy.odeint, and everything is working well, but there is something that I would like to do that I'm not sure if odeint can currently handle. The function representing the right hand side of the ODE's I generate are of the standard form that scipy needs, namely: def f(t, x, params): .... return dxdt In the ... there are many (sometimes thousands or tens of thousands of lines) intermediate variables that are used for the final calculation of dxdt. The calculations generally involve the basic arithmetic operators (+-*/**) and some trig functions (sin, cos, tan). There are several other functions of t, x, and params that I would like to evaluate *only at the integration step times specified by the user*, and I would like to be able make use of the fact that the computation of the intermediate variables has already been done for the time step that is to be returned in the x array. Here is an example of what I am trying to acheive: def f(t, x, params): z1 = ... z2 = ... (z1, z2, z3 would be in reality be more like z1, ..., z1000, all functions of t, x, params) ... zn = ... dxdt = [z1+z2**2, z2*z3-z1**3, ...] # The following only to be evaluated at time points in the user specified vector, e.g., t = arange(t_i, t_f, t_s) eqn_1 = f_1(z1,...., zn) eqn_2 = f_2(z1, ..., zn) .... eqn_n = f_n(z1, ..., zn) return dxdt the other quantities eqn1, eqn2... *after* numerical integration. But for certain applications, such as realtime control applications, it is useful to avoid repeated calculations (as would happen if I just called a function that took that state trajectory and calculated the eqn1, eqn2). I can see a couple of approaches that might work: 1) Make the time vector a global variable, and use an if statement inside of f(t, x, params) like: if t in t_span: eqn_1 = ... eqn_2 = ... .... eqn_n = ... accumulate_eqns(eqn1, ..., eqn_n) # accumulate_eqns could be either a function passed to f through pararms, or a global function 2) Make eqn_1, ..., eqn_n globals, and only integrate one time step at a time so that they can be stored at the end of each integration step. In this way, odeint would be called in a for loop that would use a time vector that had only two elements (their difference would be the step time specified by the user) I'm not 100% sure the above methods would work because I don't know if globals in the namespace which calls odeint would be viewable to globals declared in the function f. What would be really awesome would be if odeint could hold on to any variables in the f namespace when it reaches a time step that is requested, and then call another function to evaluate the output equations, but that function would have access to the state of all the variables in the f namespace *at the time point* where the state x is to be evaluated. Ideally, the extra equations wouldn't be computed except at the time points which are specified by the time vector specified by the user. Examples of use would be computing the linear/angular momentum of a system of rigid bodies and particles, center of mass of the bodies, total energy (kinetic and potential), and the value of constraints / conserved quantities during the numerical integration of the equations of motion. The equations of motion can be integrated without computing these quantities, but they are very useful for checking that the numerical integration is well behaved. For complicated systems such as spacecraft, these equations can be extremely long, and it can be very computationally advantageous to make use of these intermediate quantities rather than having to be recalculated after the numerical integration is complete. Any ideas? Sincerely, ~Luke