Consistently initialize before computing the linearization

Co-authored-by: Fredrik Bagge Carlson <baggepinnen@gmail.com>

Author: YingboMa

Project.toml

@@ -31,6 +31,7 @@ LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 MLStyle = "d8e11817-5142-5d16-987a-aa16d5891078"
 MacroTools = "1914dd2f-81c6-5fcd-8719-6d5c9610ff09"
 NaNMath = "77ba4419-2d1f-58cd-9bb1-8ffee604a2e3"
+OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
 RecursiveArrayTools = "731186ca-8d62-57ce-b412-fbd966d074cd"
 Reexport = "189a3867-3050-52da-a836-e630ba90ab69"
 RuntimeGeneratedFunctions = "7e49a35a-f44a-4d26-94aa-eba1b4ca6b47"
@@ -78,6 +79,7 @@ Latexify = "0.11, 0.12, 0.13, 0.14, 0.15, 0.16"
 MLStyle = "0.4.17"
 MacroTools = "0.5"
 NaNMath = "0.3, 1"
+OrdinaryDiffEq = "6"
 RecursiveArrayTools = "2.3"
 Reexport = "0.2, 1"
 RuntimeGeneratedFunctions = "0.5.9"
@@ -107,7 +109,6 @@ NonlinearSolve = "8913a72c-1f9b-4ce2-8d82-65094dcecaec"
 Optimization = "7f7a1694-90dd-40f0-9382-eb1efda571ba"
 OptimizationMOI = "fd9f6733-72f4-499f-8506-86b2bdd0dea1"
 OptimizationOptimJL = "36348300-93cb-4f02-beb5-3c3902f8871e"
-OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 ReferenceTests = "324d217c-45ce-50fc-942e-d289b448e8cf"
 SafeTestsets = "1bc83da4-3b8d-516f-aca4-4fe02f6d838f"

src/systems/abstractsystem.jl

@@ -1,3 +1,4 @@
+using OrdinaryDiffEq
 const SYSTEM_COUNT = Threads.Atomic{UInt}(0)
 
 get_component_type(x::AbstractSystem) = get_gui_metadata(x).type
@@ -1294,6 +1295,9 @@ function linearization_function(sys::AbstractSystem, inputs,
             if u !== nothing # Handle systems without states
                 length(sts) == length(u) ||
                     error("Number of state variables ($(length(sts))) does not match the number of input states ($(length(u)))")
+                prob = ODEProblem(fun, u, (t, t + 1), p)
+                integ = init(prob, Rodas5P())
+                u = integ.u
                 uf = SciMLBase.UJacobianWrapper(fun, t, p)
                 fg_xz = ForwardDiff.jacobian(uf, u)
                 h_xz = ForwardDiff.jacobian(let p = p, t = t