Hi all,

Just been playing around with the DifferentialEquation and ModelingToolkit systems within Turing and spotted the behaviour below.

Setting up the system and generating synthetic data:

```
using ModelingToolkit, DifferentialEquations
@parameters t
@variables y1(t) y2(t)
@parameters p1 p2
D = Differential(t)
# Define the system of equations
eqs = [
D(y1) ~ y2,
D(y2) ~ -p1 * y2 - p2 * y1
]
# Create an ODESystem
@named sys = ODESystem(eqs, t, [y1, y2], [p1, p2])
# Initial conditions
u0 = [y1 => 1.0, y2 => 0.0]
# Parameter values
p = [p1 => 1.0, p2 => 2.0]
# Time span for the solution
tspan = (0.0, 10.0)
# Convert to standard ODEProblem
sys_simplified = structural_simplify(sys)
prob = ODEProblem(sys_simplified, u0, tspan, p)
# Solve the problem
sol_orig = solve(prob)
# Synthetic data generation (for demonstration)
sim_data = sol_orig[1,:] + 0.1 * randn(length(sol_orig.t))
```

I then integrate this into Turing:

```
using Turing
# Define the model for Turing
@model function fit_model(data, prob)
# Priors for the parameters
p1 ~ Normal(1.0, 0.5)
p2 ~ Normal(2.0, 0.5)
# Solve the ODE system with the current parameter values
prob = remake(prob, p=[p1, p2]) # p = [p2, p1] ?
predicted = solve(prob, saveat=0.1)
# Likelihood function
for i in 1:length(data)
data[i] ~ Normal(predicted[i][1], 0.1)
end
end
# Run the model
model = fit_model(sim_data, prob)
chain_prep = sample(model, NUTS(), MCMCThreads(), 5, 2)
chain = sample(model, NUTS(), MCMCThreads(), 2000, 2)
```

This gives the `p1`

and `p2`

as `4.04`

and `2.018`

respectively, and lead to the following solutions, which is clearly different to that in the original set

```
# mean(chain_2[:p1])
# Initial conditions
u0 = [y1 => 1.0, y2 => 0.0]
# Parameter values
p_ch1 = [p1 => mean(chain[:p1]), p2 => mean(chain[:p2])]
p = p_ch1
# Time span for the solution
tspan = (0.0, 10.0)
# Convert to standard ODEProblem
sys_simplified = structural_simplify(sys)
prob = ODEProblem(sys_simplified, u0, tspan, p)
# Solve the problem
sol = solve(prob)
Plots.plot(sol[1, :], title = "chain_1", xlimits = ((0, 25)))
Plots.plot!(sol[2, :])
Plots.plot!(sol_orig[1, :], label = "orig y1")
Plots.plot!(sol_orig[2, :], label = "orig y2")
```

However, if I swap `p = [p1, p2]`

with `p = [p2, p1]`

in the Turing code block above, then the solution is much closer to what it should have been.

Am I missing something here? Why would I need to swap the positions of `p1`

and `p2`

? Iâ€™m using `ModelingToolkit v9.15.0`