`MonteCarloMeasurements` for a `NonlinearProblem` about modelingtoolkit.jl HOT 2 OPEN

using ModelingToolkit
using MonteCarloMeasurements

pars = @parameters a, b
vars = @variables x, y

eqs = [
    x ~ a
    y ~ x + b
]

@mtkbuild model = NonlinearSystem(eqs, vars, pars) #split=false

u0 = Dict(
    x => 1, 
    y => 2
)


function solve_wrapper(a0, b0)
    p = Dict(
        model.a => a0,
        model.b => b0
        )

    prob = NonlinearProblem(model, u0, p)
    sol = solve(prob)
    [sol[model.x], sol[model.y]]
end


a0 = 1
b0 = 1
solve_wrapper(a0, b0)

a = 1 ± 0.1
b = 1 ± 0.1
xy = @bymap solve_wrapper(a, b)

julia> xy = @bymap solve_wrapper(a, b)
2-element Vector{Particles{Float64, 2000}}:
 1.0 ± 0.1
 2.0 ± 0.14