diff --git a/Project.toml b/Project.toml
index d3ab5e4..6ff0a08 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,19 +1,23 @@
 name = "EOptInterface"
 uuid = "96d9f51d-8c41-4e2d-b1a3-7b48dea3ddd1"
 authors = ["Joseph Choi <jsphchoi@mit.edu>", "Dimitri Alston <dimitri.alston@uconn.edu>", "Pengfei Xu <pengfei.xu@uconn.edu>", "Matthew Stuber <matthew.stuber@uconn.edu>"]
-version = "0.1.1"
+version = "0.2.0"
 
 [deps]
 JuMP = "4076af6c-e467-56ae-b986-b466b2749572"
 ModelingToolkit = "961ee093-0014-501f-94e3-6117800e7a78"
+Reexport = "189a3867-3050-52da-a836-e630ba90ab69"
+SciCompDSL = "91a8cdf1-4ca6-467b-a780-87fda3fff15e"
 SymbolicUtils = "d1185830-fcd6-423d-90d6-eec64667417b"
 Symbolics = "0c5d862f-8b57-4792-8d23-62f2024744c7"
 
 [compat]
 JuMP = "1"
-ModelingToolkit = "10"
-SymbolicUtils = "3"
-Symbolics = "6"
+ModelingToolkit = "11"
+Reexport = "1"
+SciCompDSL = "1"
+SymbolicUtils = "4"
+Symbolics = "7"
 julia = "1.10"
 
 [extras]
diff --git a/docs/src/news.md b/docs/src/news.md
index 542069c..cb4961a 100644
--- a/docs/src/news.md
+++ b/docs/src/news.md
@@ -1,5 +1,9 @@
 # News for EOptInterface Releases
 
+## [v0.2.0](https://github.com/PSORLab/EOptInterface.jl/releases/tag/v0.2.0) (June 19, 2026)
+
+- Bumped ModelingToolkit.jl compatibility to v11.
+
 ## [v0.1.1](https://github.com/PSORLab/EOptInterface.jl/releases/tag/v0.1.1) (April 22, 2026)
 
 - Added notebook examples.
diff --git a/src/EOptInterface.jl b/src/EOptInterface.jl
index b6d60bf..3669e88 100644
--- a/src/EOptInterface.jl
+++ b/src/EOptInterface.jl
@@ -16,9 +16,13 @@ module EOptInterface
     import JuMP
     import ModelingToolkit
     import ModelingToolkit: t_nounits as t, D_nounits as D
+    import Reexport
     import Symbolics
     import SymbolicUtils
 
+    # Reexports
+    Reexport.@reexport using SciCompDSL
+
     # Include base functions
     include("basefuncs.jl")
 
diff --git a/src/basefuncs.jl b/src/basefuncs.jl
index 93aad10..8865532 100644
--- a/src/basefuncs.jl
+++ b/src/basefuncs.jl
@@ -12,7 +12,7 @@
 ################################################################################
 
 function mtk_generate_model_equations(sys::ModelingToolkit.System)
-    param_dict = copy(ModelingToolkit.defaults(sys))
+    param_dict = copy(ModelingToolkit.initial_conditions(sys))
     h = []
     for i in eachindex(ModelingToolkit.unknowns(sys))
         expr = (ModelingToolkit.full_equations(ModelingToolkit.expand_connections(sys))[i].rhs 
@@ -31,7 +31,7 @@ function mtk_generate_model_equations(sys::ModelingToolkit.System)
 end
 
 function mtk_generate_reduced_expression(expr::Symbolics.Num, sys::ModelingToolkit.System)
-    sub_dict = ModelingToolkit.defaults(sys)
+    sub_dict = ModelingToolkit.initial_conditions(sys)
     for eqn in ModelingToolkit.observed(sys)
         sub_dict[eqn.lhs] = eqn.rhs
     end
diff --git a/src/userfuncs.jl b/src/userfuncs.jl
index 9213779..26497ae 100644
--- a/src/userfuncs.jl
+++ b/src/userfuncs.jl
@@ -20,7 +20,7 @@ Returns the decision variables for an optimization problem formulated from a Mod
 function decision_vars(sys::ModelingToolkit.System)
     return [
         ModelingToolkit.unknowns(sys); 
-        setdiff(ModelingToolkit.parameters(sys), keys(ModelingToolkit.defaults(sys)))
+        setdiff(ModelingToolkit.parameters(sys), keys(ModelingToolkit.initial_conditions(sys)))
         ]
 end
 
@@ -70,13 +70,15 @@ function register_odesystem(model::JuMP.Model, odesys::ModelingToolkit.System, t
     N = Int(floor((tspan[2] - tspan[1])/tstep)) + 1
     # Number of ODE variables
     V = length(ModelingToolkit.unknowns(odesys))
-    param_dict = copy(ModelingToolkit.defaults(odesys))
+    param_dict = copy(ModelingToolkit.initial_conditions(odesys).dict)
     for var in ModelingToolkit.unknowns(odesys)
         pop!(param_dict, var)
     end
     dx = []
     for j in 1:V
         dxj_expr = ModelingToolkit.full_equations(odesys)[j].rhs
+        dxj_expr = SymbolicUtils.substitute(dxj_expr, ModelingToolkit.bindings(odesys))
+        dxj_expr = SymbolicUtils.substitute(dxj_expr, ModelingToolkit.bindings(odesys))
         # Fully substitute parameters with default values
         while ~isempty(intersect(Symbolics.get_variables(dxj_expr), keys(param_dict)))
             dxj_expr = SymbolicUtils.substitute(dxj_expr, param_dict)
@@ -91,7 +93,7 @@ function register_odesystem(model::JuMP.Model, odesys::ModelingToolkit.System, t
     ps = JuMP.all_variables(model)[end-length(setdiff(EOptInterface.decision_vars(odesys), ModelingToolkit.unknowns(odesys)))+1:end]
     xs = reshape(setdiff(JuMP.all_variables(model), ps), V, N)
     # Extract initial conditions from the ModelingToolkit system and fix them in the JuMP model for x[1:V,1]
-    JuMP.fix.(xs[:,1], [ModelingToolkit.defaults(odesys)[ModelingToolkit.unknowns(odesys)[i]] for i in eachindex(ModelingToolkit.unknowns(odesys))], force=true)
+    JuMP.fix.(xs[:,1], [ModelingToolkit.initial_conditions(odesys)[ModelingToolkit.unknowns(odesys)[i]].val for i in eachindex(ModelingToolkit.unknowns(odesys))], force=true)
     # Formulate JuMP constraints based on chosen ODE discretization method
     if integrator == "EE"
         JuMP.@constraint(model, [j in 1:V, i in 1:(N-1)], xs[j,i+1] == xs[j,i] + tstep*dx[j](xs[:,i]..., ps...))
@@ -108,7 +110,7 @@ Returns a dictionary of optimal solution values for the observed variables of an
 """
 function full_solution(model::JuMP.Model, sys::ModelingToolkit.System)
     vars = EOptInterface.decision_vars(sys)
-    sub_dict = ModelingToolkit.defaults(sys)
+    sub_dict = ModelingToolkit.initial_conditions(sys)
     for i in eachindex(vars)
         sub_dict[vars[i]] = JuMP.value.(JuMP.all_variables(model)[i])
     end
diff --git a/test/runtests.jl b/test/runtests.jl
index db46732..fef9a09 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -3,6 +3,7 @@ using Ipopt
 using JuMP
 using ModelingToolkit
 using ModelingToolkit: t_nounits as t, D_nounits as D
+using SciCompDSL
 using Test
 
 @testset "Algebraic Model" begin