qojulia · oameye · Sep 30, 2025 · Jun 6, 2025 · Jun 21, 2025 · Jun 21, 2025
diff --git a/.github/workflows/Documenter.yml b/.github/workflows/Documenter.yml
@@ -40,7 +40,7 @@ jobs:
       - uses: actions/checkout@v4
       - uses: julia-actions/setup-julia@v2
         with:
-          version: 'lts'
+          version: '1'
       - uses: julia-actions/cache@v2
       - uses: julia-actions/julia-buildpkg@v1
       - uses: julia-actions/julia-docdeploy@v1

diff --git a/.github/workflows/Tests.yml b/.github/workflows/Tests.yml
@@ -46,7 +46,7 @@ jobs:
       fail-fast: false
       matrix:
         version:
-          - 'pre'
+          # - 'pre'
           - 'lts'
           - '1'
         os:

diff --git a/Project.toml b/Project.toml
@@ -27,14 +27,14 @@ LaTeXStrings = "1"
 Latexify = "0.13, 0.14, 0.15, 0.16"
 LinearAlgebra = "1.10"
 MacroTools = "0.5"
-ModelingToolkit = "~9.70"
+ModelingToolkit = "10"
 OrdinaryDiffEq = "6"
 OrderedCollections = "1.8"
 QuantumOptics = "1"
 QuantumOpticsBase = "0.4, 0.5"
 Random = "1.10"
 SciMLBase = "1, 2"
-SecondQuantizedAlgebra = "0.3.1"
+SecondQuantizedAlgebra = "=0.4.1"
 SteadyStateDiffEq = "2"
 SymbolicUtils = "3.6"
 Symbolics = "6"

diff --git a/README.md b/README.md
@@ -44,10 +44,10 @@ rates = [κ,γ,ν]
 eqs = meanfield([a,σ(:g,:e),σ(:e,:e)], H, J; rates=rates, order=1)
 
 using ModelingToolkit, OrdinaryDiffEq
-@named sys = ODESystem(eqs)
-p0 = (Δ=>0, g=>1.5, κ=>1, γ=>0.25, ν=>4)
-u0 = ComplexF64[1e-2, 0, 0]
-prob = ODEProblem(sys,u0,(0.0,50.0),p0)
+@named sys = System(eqs)
+p0 = Dict(Δ=>0, g=>1.5, κ=>1, γ=>0.25, ν=>4)
+u0 = Dict(unknowns(sys) => ComplexF64[1e-2, 0, 0])
+prob = ODEProblem(sys,merge(u0, p0),(0.0,50.0))
 sol = solve(prob,RK4())
 
 using Plots
@@ -58,7 +58,7 @@ plot(sol.t, n, xlabel="t", label="n")
 ![photon-number](https://user-images.githubusercontent.com/18166442/114183684-3ae76080-9944-11eb-9d21-94bf4069bb60.png)
 
 
-The above code implements the Jaynes-Cummings Hamiltonian describing an optical cavity mode that couples to a two-level atom. Additionally, the decay processes are specified. Then, mean-field equations for the average values of the operators `[a,σ(:g,:e),σ(:e,:e)]` are derived and expanded to first order (average values of products are factorized). For the numerical solution an `ODESystem` (from [ModelingToolkit.jl](https://github.com/SciML/ModelingToolkit.jl)) is created and solved with the [OrdinaryDiffEq.jl](https://github.com/SciML/OrdinaryDiffEq.jl) library. Finally, the time dynamics of the photon number `n` is plotted.
+The above code implements the Jaynes-Cummings Hamiltonian describing an optical cavity mode that couples to a two-level atom. Additionally, the decay processes are specified. Then, mean-field equations for the average values of the operators `[a,σ(:g,:e),σ(:e,:e)]` are derived and expanded to first order (average values of products are factorized). For the numerical solution an `System` (from [ModelingToolkit.jl](https://github.com/SciML/ModelingToolkit.jl)) is created and solved with the [OrdinaryDiffEq.jl](https://github.com/SciML/OrdinaryDiffEq.jl) library. Finally, the time dynamics of the photon number `n` is plotted.
 
 
 ## Citing

diff --git a/docs/Project.toml b/docs/Project.toml
@@ -17,11 +17,10 @@ Symbolics = "0c5d862f-8b57-4792-8d23-62f2024744c7"
 Documenter = "1"
 Latexify = "0.13, 0.14, 0.15, 0.16"
 MacroTools = "0.5"
-ModelingToolkit = "9"
 OrdinaryDiffEq = "6"
 Plots = "1"
 QuantumOptics = "1"
 QuantumOpticsBase = "0.4.21"
 StochasticDiffEq = "6"
-SymbolicUtils = "3.6.0 - 3.20.0"
+SymbolicUtils = "3.26.0"
 Symbolics = "6"
diff --git a/docs/src/correlation.md b/docs/src/correlation.md
@@ -47,27 +47,27 @@ nothing # hide
 ```
 When the [`CorrelationFunction`](@ref) is constructed, an additional Hilbert space is added internally which represents the system at the time ``t``. In our case, this means that another [`FockSpace`](@ref) is added. Note that all operators involved in the correlation function are defined on the [`ProductSpace`](@ref) including this additional Hilbert space.
 
-The equation for ``g(t,\tau)`` is now stored in the first entry of `c.de`. To solve the above numerically, we need to convert to an `ODESystem` and solve numerically.
+The equation for ``g(t,\tau)`` is now stored in the first entry of `c.de`. To solve the above numerically, we need to convert to an `System` and solve numerically.
 ```@example correlation
 using ModelingToolkit, OrdinaryDiffEq
 
-@named sys = ODESystem(me)
+@named sys = System(me)
 n0 = 20.0 # Initial number of photons in the cavity
-u0 = [n0]
-p0 = (ωc => 1, κ => 1)
-prob = ODEProblem(sys,u0,(0.0,2.0),p0) # End time not in steady state
+n_avg = unknowns(sys) |> first
+p0 = (ωc => 1, κ => 1, n_avg => n0) # Initial values and parameters
+prob = ODEProblem(sys,Dict(p0),(0.0,2.0)) # End time not in steady state
 sol = solve(prob,RK4())
 nothing # hide
 ```
 Numerically computing the correlation function works in the same way. Note, the initial state of the correlation function depends on the final state of the system. However, in general it does not depend on *all* the final values of the system. The correct values can be picked out automatically using the [`correlation_u0`](@ref) function.
 ```@example correlation
-@named csys = ODESystem(c)
+@named csys = System(c)
 u0_c = correlation_u0(c, sol.u[end])
-prob_c = ODEProblem(csys,u0_c,(0.0,10.0),p0)
+prob_c = ODEProblem(csys,merge(Dict(u0_c),Dict(p0)),(0.0,10.0))
 sol_c = solve(prob_c,RK4(),save_idxs=1)
 nothing # hide
 ```
-Finally, lets check our numerical solution against the analytic one obtained above:
+Finally, let's check our numerical solution against the analytic one obtained above:
 ```@example correlation
 using Test # hide
 g_analytic(τ) = @. sol.u[end] * exp((im*p0[1][2]-0.5p0[2][2])*τ)

diff --git a/docs/src/examples/cavity_antiresonance_indexed.md b/docs/src/examples/cavity_antiresonance_indexed.md
@@ -101,7 +101,7 @@ To create the equations for a specific number of atoms we use the function $\tex
 ```@example antiresonance_indexed
 N_ = 2
 eqs_ = evaluate(eqs;limits=(N=>N_))
-@named sys = ODESystem(eqs_)
+@named sys = System(eqs_)
 nothing # hide
 ```
 
@@ -149,7 +149,8 @@ for i=1:length(Δ_ls)
     p0_ = [Δc_i; η_; Δa_i; κ_; gi_; Γij_; Ωij_]
 
     # create (remake) new ODEProblem
-    prob_  = ODEProblem(sys,u0,(0.0, 20),ps.=>p0_)
+    dict = merge(Dict(unknowns(sys) .=> u0), Dict(ps .=> p0_))
+    prob_  = ODEProblem(sys,dict,(0.0, 20))
     sol_ = solve(prob_, Tsit5())
     n_ls[i] = abs2(sol_[a][end])
 end

diff --git a/docs/src/examples/excitation-transport-chain.md b/docs/src/examples/excitation-transport-chain.md
@@ -64,12 +64,12 @@ Rather, we only need to derive the equations once and substitute the noisy posit
 eqs = meanfield(σ(:g,:e,1),H,c_ops;rates=[γ for i=1:N],order=2)
 complete!(eqs)  # complete the set
 
-# Generate the ODESystem
-@named sys = ODESystem(eqs)
+# Generate the System
+@named sys = System(eqs)
 nothing # hide
 ```
 
-Once we have our set of equations and converted it to an `ODESystem` we are ready to solve for the dynamics.
+Once we have our set of equations and converted it to an `System` we are ready to solve for the dynamics.
 First, let's have a look at the excitation transport for perfectly positioned atoms.
 We assume an equidistant chain, were neighboring atoms are separated by a distance $d$.
 
@@ -81,7 +81,8 @@ p = [γ => 1.0; Δ => 0.0; Ω => 2.0; J0 => 1.25; x .=> x0;]
 
 # Create ODEProblem
 u0 = zeros(ComplexF64, length(eqs))  # initial state -- all atoms in the ground state
-prob = ODEProblem(sys,u0,(0.0,15.0),p)
+dict = merge(Dict(unknowns(sys) .=> u0), Dict(p))
+prob = ODEProblem(sys,dict,(0.0,15.0))
 
 # Solve
 sol = solve(prob,RK4())
@@ -112,7 +113,8 @@ function prob_func(prob,i,repeat)
     x_ = x0 .+ s.*randn(N)
     p_ = [γ => 1.0; Δ => 0.0; Ω => 2.0; J0 => 1.25; x .=> x_;]
     # Return new ODEProblem
-    return ODEProblem(sys,u0,(0.0,15.0),p_)
+    dict = merge(Dict(unknowns(sys) .=> u0), Dict(p_))
+    return ODEProblem(sys,dict,(0.0,15.0))
 end
 
 trajectories = 20

diff --git a/docs/src/examples/filter-cavity_indexed.md b/docs/src/examples/filter-cavity_indexed.md
@@ -99,7 +99,7 @@ To calculate the dynamics of the system we create a system of ordinary different
 
 
 ```@example filter_cavity_indexed
-@named sys = ODESystem(eqs_eval)
+@named sys = System(eqs_eval)
 nothing # hide
 ```
 
@@ -124,7 +124,8 @@ gf_ = 0.1Γ_
 ps = [Γ, κ, g, κf, gf, R, [δ(i) for i=1:M_]..., Δ, ν, N]
 p0 = [Γ_, κ_, g_, κf_, gf_, R_, δ_ls..., Δ_, ν_, N_]
 
-prob = ODEProblem(sys,u0,(0.0, 10.0/κf_), ps.=>p0)
+dict = merge(Dict(unknowns(sys) .=> u0), Dict(ps .=> p0))
+prob = ODEProblem(sys,dict,(0.0, 10.0/κf_))
 nothing # hide
 ```
 

diff --git a/docs/src/examples/heterodyne_detection.md b/docs/src/examples/heterodyne_detection.md
@@ -88,7 +88,7 @@ scaled_eqs = scale(eqs_c)
 ```
 
 # Deterministic time evolution
-Here we define the actual values for the system parameters. We then show that the deterministic time evolution without the noise terms is still accessible by using the constructor ODESystem for the stochastic system of equations and the syntax for the simulation of the time evolution is as usual.
+Here we define the actual values for the system parameters. We then show that the deterministic time evolution without the noise terms is still accessible by using the constructor System for the stochastic system of equations and the syntax for the simulation of the time evolution is as usual.
 
 
 ```julia
@@ -98,9 +98,10 @@ g_ = 6.531*10^3; ωl_ = 2.0 * π * 10^3; t0=0.0;t1=20e-6
 p = [N,ωa,γ,η,χ,ωc,κ,g,ξ,ωl]
 p0 = [N_,ωa_,γ_,η_,χ_,ωc_,κ_,g_,ξ_,ωl_]
 
-@named sys = ODESystem(scaled_eqs)
+@named sys = System(scaled_eqs)
 u0 = zeros(ComplexF64,length(scaled_eqs.equations))
-prob = ODEProblem(sys, u0,(0.0,1e-3),p.=>p0)
+dict = merge(Dict(unknowns(sys) .=> u0), Dict(p .=> p0))
+prob = ODEProblem(sys, dict,(0.0,1e-3))
 sol_det = solve(prob,RK4(),dt=1e-9)
 
 plot(sol_det.t .* 1000, map(x -> real(x[2]), sol_det.u))
@@ -110,7 +111,7 @@ xlabel("Time (ms)")
 ```
 
 # Stochastic time evolution
-The stochastic time evolution is accessible via the constructor SDESystem, whose syntax is exactly the same as for the ODESystem, but with keyword args as defined in https://docs.sciml.ai/DiffEqDocs/stable/tutorials/sde_example/. We then need to provide a noise process for the measurement. If the noise is white the appropriate noise process is a Wiener process. The SDEProblem is then constructed just as the ODEProblem, but with an additional noise argument.
+The stochastic time evolution is accessible via the constructor SDESystem, whose syntax is exactly the same as for the System, but with keyword args as defined in https://docs.sciml.ai/DiffEqDocs/stable/tutorials/sde_example/. We then need to provide a noise process for the measurement. If the noise is white the appropriate noise process is a Wiener process. The SDEProblem is then constructed just as the ODEProblem, but with an additional noise argument.
 
 We can then make use of the EnsembleProblem which automatically runs multiple instances of the stochastic equations of motion. The number of trajectories can then be set in the solve call. See the documentation cited above for more details of the function calls here.
 

diff --git a/docs/src/examples/jupyter_notebooks/cavity_antiresonance.ipynb b/docs/src/examples/jupyter_notebooks/cavity_antiresonance.ipynb
@@ -157,7 +157,7 @@
    "outputs": [],
    "source": [
     "complete!(eqs)\n",
-    "@named sys = ODESystem(eqs)\n",
+    "@named sys = System(eqs)\n",
     "nothing # hide"
    ]
   },

diff --git a/docs/src/examples/jupyter_notebooks/cavity_antiresonance_indexed.ipynb b/docs/src/examples/jupyter_notebooks/cavity_antiresonance_indexed.ipynb
@@ -186,13 +186,13 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 6,
+   "execution_count": null,
    "metadata": {},
    "outputs": [],
    "source": [
     "N_ = 2\n",
     "eqs_ = evaluate(eqs;limits=(N=>N_))\n",
-    "@named sys = ODESystem(eqs_)\n",
+    "@named sys = System(eqs_)\n",
     "nothing # hide"
    ]
   },

diff --git a/docs/src/examples/jupyter_notebooks/filter-cavity_indexed.ipynb b/docs/src/examples/jupyter_notebooks/filter-cavity_indexed.ipynb
@@ -201,7 +201,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "@named sys = ODESystem(eqs_eval)\n",
+    "@named sys = System(eqs_eval)\n",
     "nothing # hide"
    ]
   },

diff --git a/docs/src/examples/jupyter_notebooks/heterodyne_detection.ipynb b/docs/src/examples/jupyter_notebooks/heterodyne_detection.ipynb
@@ -395,7 +395,7 @@
    "metadata": {},
    "source": [
     "# Deterministic time evolution\n",
-    "Here we define the actual values for the system parameters. We then show that the deterministic time evolution without the noise terms is still accessible by using the constructor ODESystem for the stochastic system of equations and the syntax for the simulation of the time evolution is as usual."
+    "Here we define the actual values for the system parameters. We then show that the deterministic time evolution without the noise terms is still accessible by using the constructor System for the stochastic system of equations and the syntax for the simulation of the time evolution is as usual."
    ]
   },
   {
@@ -423,7 +423,7 @@
     "p = [N,ωa,γ,η,χ,ωc,κ,g,ξ,ωl]\n",
     "p0 = [N_,ωa_,γ_,η_,χ_,ωc_,κ_,g_,ξ_,ωl_]\n",
     "\n",
-    "@named sys = ODESystem(scaled_eqs)\n",
+    "@named sys = System(scaled_eqs)\n",
     "u0 = zeros(ComplexF64,length(scaled_eqs.equations))\n",
     "prob = ODEProblem(sys, u0,(0.0,1e-3),p.=>p0)\n",
     "sol_det = solve(prob,RK4(),dt=1e-9)\n",
@@ -439,7 +439,7 @@
    "metadata": {},
    "source": [
     "# Stochastic time evolution\n",
-    "The stochastic time evolution is accessible via the constructor SDESystem, whose syntax is exactly the same as for the ODESystem, but with keyword args as defined in https://docs.sciml.ai/DiffEqDocs/stable/tutorials/sde_example/. We then need to provide a noise process for the measurement. If the noise is white the appropriate noise process is a Wiener process. The SDEProblem is then constructed just as the ODEProblem, but with an additional noise argument.\n",
+    "The stochastic time evolution is accessible via the constructor SDESystem, whose syntax is exactly the same as for the System, but with keyword args as defined in https://docs.sciml.ai/DiffEqDocs/stable/tutorials/sde_example/. We then need to provide a noise process for the measurement. If the noise is white the appropriate noise process is a Wiener process. The SDEProblem is then constructed just as the ODEProblem, but with an additional noise argument.\n",
     "\n",
     "We can then make use of the EnsembleProblem which automatically runs multiple instances of the stochastic equations of motion. The number of trajectories can then be set in the solve call. See the documentation cited above for more details of the function calls here."
    ]

diff --git a/docs/src/examples/jupyter_notebooks/many_atom_laser.ipynb b/docs/src/examples/jupyter_notebooks/many_atom_laser.ipynb
@@ -143,8 +143,8 @@
     "To calculate the time evolution we create a Julia function which can be used by DifferentialEquations.jl to solve the set of ordinary differential equations.\n",
     "\n",
     "```@example 3-level-laser\n",
-    "# Build an ODESystem out of the MeanfieldEquations\n",
-    "@named sys = ODESystem(me_comp)\n",
+    "# Build an System out of the MeanfieldEquations\n",
+    "@named sys = System(me_comp)\n",
     "nothing # hide"
    ]
   },

diff --git a/docs/src/examples/jupyter_notebooks/mollow triplet.ipynb b/docs/src/examples/jupyter_notebooks/mollow triplet.ipynb
@@ -182,7 +182,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "@named sys = ODESystem(eqs)\n",
+    "@named sys = System(eqs)\n",
     "\n",
     "p0 = (0.0,2.0,1.0)\n",
     "u0 = zeros(ComplexF64, 2)\n",

diff --git a/docs/src/examples/jupyter_notebooks/noisy excitation transport.ipynb b/docs/src/examples/jupyter_notebooks/noisy excitation transport.ipynb
@@ -96,16 +96,16 @@
     "eqs = meanfield(σ(:g,:e,1),H,c_ops;rates=[γ for i=1:N],order=2)\n",
     "complete!(eqs)  # complete the set\n",
     "\n",
-    "# Generate the ODESystem\n",
-    "@named sys = ODESystem(eqs)\n",
+    "# Generate the System\n",
+    "@named sys = System(eqs)\n",
     "nothing #hide"
    ]
   },
   {
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "Once we have our set of equations and converted it to an `ODESystem` we are ready to solve for the dynamics.\n",
+    "Once we have our set of equations and converted it to an `System` we are ready to solve for the dynamics.\n",
     "First, let's have a look at the excitation transport for perfectly positioned atoms.\n",
     "We assume an equidistant chain, were neighboring atoms are separated by a distance $d$."
    ]

diff --git a/docs/src/examples/jupyter_notebooks/optomechanical cooling.ipynb b/docs/src/examples/jupyter_notebooks/optomechanical cooling.ipynb
@@ -117,7 +117,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "@named sys = ODESystem(eqs_completed)\n",
+    "@named sys = System(eqs_completed)\n",
     "nothing #hide"
    ]
   },

diff --git a/docs/src/examples/jupyter_notebooks/ramsey spectroscopy.ipynb b/docs/src/examples/jupyter_notebooks/ramsey spectroscopy.ipynb
@@ -103,7 +103,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "@named sys = ODESystem(eqs)\n",
+    "@named sys = System(eqs)\n",
     "\n",
     "# Parameter\n",
     "Γ_ = 1.0\n",

diff --git a/docs/src/examples/jupyter_notebooks/spectrum of single atom laser.ipynb b/docs/src/examples/jupyter_notebooks/spectrum of single atom laser.ipynb
@@ -173,7 +173,7 @@
    "source": [
     "# Numerical solution\n",
     "ps = (Δ, g, γ, κ, ν)\n",
-    "@named sys = ODESystem(eqs)\n",
+    "@named sys = System(eqs)\n",
     "u0 = zeros(ComplexF64, length(eqs))\n",
     "p0 = (1.0, 1.5, 0.25, 1, 4)\n",
     "prob = ODEProblem(sys,u0,(0.0,10.0),ps.=>p0)\n",
@@ -196,7 +196,7 @@
    "outputs": [],
    "source": [
     "# Time evolution of correlation function\n",
-    "@named csys = ODESystem(c)\n",
+    "@named csys = System(c)\n",
     "u0_c = correlation_u0(c,sol.u[end])\n",
     "p0_c = correlation_p0(c,sol.u[end],ps.=>p0)\n",
     "prob_c = ODEProblem(csys,u0_c,(0.0,500.0),p0_c)\n",

diff --git a/docs/src/examples/jupyter_notebooks/superradiant_laser.ipynb b/docs/src/examples/jupyter_notebooks/superradiant_laser.ipynb
@@ -182,7 +182,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "@named sys = ODESystem(eqs_c)\n",
+    "@named sys = System(eqs_c)\n",
     "nothing # hide"
    ]
   },

diff --git a/docs/src/examples/jupyter_notebooks/superradiant_laser_indexed.ipynb b/docs/src/examples/jupyter_notebooks/superradiant_laser_indexed.ipynb
@@ -256,7 +256,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "@named sys = ODESystem(eqs_sc)\n",
+    "@named sys = System(eqs_sc)\n",
     "nothing # hide"
    ]
   },

diff --git a/docs/src/examples/jupyter_notebooks/unique_squeezing.ipynb b/docs/src/examples/jupyter_notebooks/unique_squeezing.ipynb
@@ -209,7 +209,7 @@
    "outputs": [],
    "source": [
     "# symbolic ordinary differential equation system\n",
-    "@named sys = ODESystem(eqs_sc)\n",
+    "@named sys = System(eqs_sc)\n",
     "\n",
     "# initial state \n",
     "u0 = zeros(ComplexF64, length(eqs_sc));\n",
@@ -338,7 +338,7 @@
    "outputs": [],
    "source": [
     "# symbolic ordinary differential equation system\n",
-    "@named sys_a = ODESystem(eqs_a)\n",
+    "@named sys_a = System(eqs_a)\n",
     "\n",
     "# initial state \n",
     "u0_a = zeros(ComplexF64, length(eqs_a))\n",

diff --git a/docs/src/examples/jupyter_notebooks/waveguide.ipynb b/docs/src/examples/jupyter_notebooks/waveguide.ipynb
@@ -235,7 +235,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "@named sys = ODESystem(eqs)\n",
+    "@named sys = System(eqs)\n",
     "prob = ODEProblem(sys,u0,(0.0, 8e-3), ps.=>p0)\n",
     "sol = solve(prob,Tsit5(),abstol=1e-6,reltol=1e-6)\n",
     "nothing # hide"
-Original file line number
+Diff line change
@@ Expand Up / @@ -46,7 +46,7 @@ jobs: @@
           fail-fast: false
           matrix:
             version:
-              - 'pre'
+              # - 'pre'
               - 'lts'
               - '1'
             os:
@@ Expand Down @@