This function is intended to be called repeatedly until it converges to some value of `Xnew`

:

```
using DifferentialEquations
using ControlSystemsBase.MatrixEquations: lyapd
function newton_iterate(X::Union{ODESolution,Function},
mathieu::GeneralizedMathieu)
lqrparams = LQRParamsMathieu()
T = 2*pi/mathieu.ω
M = B * inv(lqrparams.r) * B'
ϕ = state_transition(mathieu, M, X)
# Find Xnew_init = X_{n+1}(0) == X_{n+1}(T)
integrand(τ) = ϕ(τ)' * (lqrparams.Q + X(τ) * M * X(τ)) * ϕ(τ)
Xnew_init = lyapd(ϕ(T)', ∫(integrand, 0, T))
# Find Xnew(t) = X_{n+1}(t) by integrating the periodic Lyapunov equation
L_operator(Xnp1, _, t) =
-(A(t, mathieu, M, X)' * Xnp1 + Xnp1 * A(t, mathieu, M, X) + (X(t) * M * X(t) + lqrparams.Q))
PLEProb = ODEProblem(L_operator, Xnew_init, (T, 0), nothing)
Xnew = solve(PLEProb, Tsit5(), saveat=T/50.0)
end
```

Unfortunately,

- it compiles on every call, and
- around the sixth iteration it consumed 10+ GB of RAM and hung indefinitely (1+ hour).

I’m trying to understand why it compiles on every call. I suspect the RAM issue is related. Here is a summary of what it does:

- Accept as input a 1-arg function or ODE solution
`X`

. -
`X`

is used as a parameter to find two intermediate quantites:`ϕ`

and`Xnew_init`

. -
`ϕ(t)`

is the solution to an ODE that depends on`X(t)`

.`Xnew_init`

is found by integrating`X(t)`

and solving an algebraic equation. - Finally,
`ϕ(t)`

,`Xnew_init`

and`X(t)`

are used in the inputs to an ODEProblem whose solution is`Xnew`

. - The function can then be called with
`Xnew`

as input, starting the process over. (It’s a Newton-Raphson iteration step applied to a time-varying Differential Ricatti Equation.)

The other quantities (`B`

, `mathieu`

, `lqrparams`

etc) are just constants defined in the function or in the calling environment. `∫`

is a simple Gaussian quadrature-based integrator. I can supply more details if required.