diff --git a/Project.toml b/Project.toml
index b4c7d9e..039a3a1 100644
--- a/Project.toml
+++ b/Project.toml
@@ -23,8 +23,10 @@ Symbolics = "7"
 julia = "1.10"
 
 [extras]
+CSV = "336ed68f-0bac-5ca0-87d4-7b16caf5d00b"
+DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"
 Ipopt = "b6b21f68-93f8-5de0-b562-5493be1d77c9"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
-test = ["Ipopt", "Test"]
+test = ["CSV", "DataFrames", "Ipopt", "Test"]
diff --git a/README.md b/README.md
index 07a64f1..010fa7b 100644
--- a/README.md
+++ b/README.md
@@ -23,9 +23,9 @@ An optimal reactor-separator-recycle process design problem originally presented
 A nonlinear kinetic parameter estimation problem originally described by [[5](#references)] is used to demonstrate the use of `register_odesystem` to formulate and solve an ODE system using [`Ipopt`](https://github.com/coin-or/ipopt) [[6](#references)].
 
 ## References
-1. Ma, Y., Gowda, S., Anantharaman, R., Laughman, C., Shah, V., and Rackauckas, C. ModelingToolkit: A Composable Graph Transformation System For Equation-Based Modeling. (2022). DOI: [10.48550/arXiv.2103.05244)](https://doi.org/10.48550/arXiv.2103.05244)
+1. Ma, Y., Gowda, S., Anantharaman, R., Laughman, C., Shah, V., and Rackauckas, C. ModelingToolkit: A Composable Graph Transformation System For Equation-Based Modeling. (2022). DOI: [10.48550/arXiv.2103.05244](https://doi.org/10.48550/arXiv.2103.05244)
 2. Lubin, M., Dowson, O., Dias Garcia, J., Huchette, J., Legat, B., and Vielma, J.P. JuMP 1.0: recent improvements to a modeling language for mathematical optimization. *Mathematical Programming Computation.* 15, 581-589 (2023). DOI: [10.1007/s12532-023-00239-3](https://doi.org/10.1007/s12532-023-00239-3)
 3. Kokossis, A.C. and Floudas, C.A. Synthesis of isothermal reactor-separator-recycle systems. *Chemical Engineering Science.* 46, 1361-1383 (1991). DOI: [10.1016/0009-2509(91)85063-4](https://doi.org/10.1016/0009-2509(91)85063-4)
-4. Wilhelm, M. E. and Stuber, M.D. EAGO.jl: easy advanced global optimization in Julia. *Optimization Methods and Software.* 37(2), 425-450 (2022). DOI: [10.1080/10556788.2020.1786566](https://doi.org/10.1080/10556788.2020.1786566)
+4. Wilhelm, M. E. and Stuber, M.D. EAGO.jl: easy advanced global optimization in Julia. *Optimization Methods & Software.* 37(2), 425-450 (2022). DOI: [10.1080/10556788.2020.1786566](https://doi.org/10.1080/10556788.2020.1786566)
 5. Taylor, J.W. Direct measurement and analysis of cyclohexadienyl oxidation. Ph.D. thesis, Massachusetts Institute of Technology. (2005). URL: http://hdl.handle.net/1721.1/33716
 6. Wächter, A. and Biegler, L.T. On the implementation of an interior-point filter line-search algorithm for large-scale nonlinear programming. *Mathematical Programming.* 106, 25-57 (2006). DOI: [10.1007/s10107-004-0559-y](https://doi.org/10.1007/s10107-004-0559-y)
diff --git a/examples/algebraic_model.ipynb b/examples/algebraic_model.ipynb
index e11537c..f534bed 100644
--- a/examples/algebraic_model.ipynb
+++ b/examples/algebraic_model.ipynb
@@ -43,12 +43,15 @@
   },
   {
    "cell_type": "code",
-   "execution_count": null,
+   "execution_count": 1,
    "id": "62f7a7b2",
    "metadata": {},
    "outputs": [],
    "source": [
-    "using ModelingToolkit, JuMP, EAGO, EOptInterface\n",
+    "using EAGO\n",
+    "using EOptInterface\n",
+    "using JuMP\n",
+    "using ModelingToolkit\n",
     "using ModelingToolkit: t_nounits as t, D_nounits as D"
    ]
   },
@@ -57,7 +60,7 @@
    "id": "f82c197c",
    "metadata": {},
    "source": [
-    "The system is modeling in `ModelingToolkit` by defining a `ModelingToolkit.@mtkmodel` for each component (or unit operation) illustrated in the figure and a `ModelingToolkit.@connector` is defined for the variables in the stream."
+    "The system is modeling in `ModelingToolkit` [[2](#references)] by defining a `ModelingToolkit.@mtkmodel` for each component (or unit operation) illustrated in the figure and a `ModelingToolkit.@connector` is defined for the variables in the stream."
    ]
   },
   {
@@ -274,7 +277,7 @@
    "id": "9777cd54",
    "metadata": {},
    "source": [
-    "We can then define the `JuMP.Model` using an appropriate solver of choice. To solve this problem to global optimality, we have chosen to use `EAGO` [[2](#references)]."
+    "We can then define the `JuMP.Model` [[3](#references)] using an appropriate solver of choice. To solve this problem to global optimality, we have chosen to use `EAGO` [[4](#references)]."
    ]
   },
   {
@@ -319,11 +322,11 @@
     {
      "data": {
       "text/plain": [
-       "6-element Vector{Any}:\n",
-       " sep1₊in₊F(t)\n",
-       " sep1₊in₊y_B(t)\n",
-       " sep1₊in₊y_C(t)\n",
-       " sep1₊outL₊y_C(t)\n",
+       "6-element Vector{SymbolicUtils.BasicSymbolicImpl.var\"typeof(BasicSymbolicImpl)\"{SymReal}}:\n",
+       " mixer₊in2₊F(t)\n",
+       " cstr₊out₊y_C(t)\n",
+       " cstr₊out₊y_B(t)\n",
+       " sep1₊outL₊y_B(t)\n",
        " influent₊F\n",
        " cstr₊V"
       ]
@@ -388,18 +391,6 @@
    "id": "1d4b16a1",
    "metadata": {},
    "outputs": [
-    {
-     "name": "stderr",
-     "output_type": "stream",
-     "text": [
-      "┌ Warning: At least one branching variable is unbounded. This will interfere with EAGO's global\n",
-      "│ optimization routine and may cause unexpected results. Bounds have been automatically\n",
-      "│ generated at +/- 1E6 for all unbounded variables, but tighter user-defined bounds are\n",
-      "│ highly recommended. To disable this warning and the automatic generation of bounds, use\n",
-      "│ the option `unbounded_check = false`.\n",
-      "└ @ EAGO C:\\Users\\ybr24001\\.julia\\packages\\EAGO\\WUXx0\\src\\eago_optimizer\\optimize\\nonconvex\\stack_management.jl:255\n"
-     ]
-    },
     {
      "name": "stdout",
      "output_type": "stream",
@@ -414,36 +405,36 @@
       "-----------------------------------------------------------------------------------------------------------------\n",
       "| Iteration # |    Nodes    | Lower Bound | Upper Bound |     Gap     |    Ratio    |    Timer    |  Time Left  |\n",
       "-----------------------------------------------------------------------------------------------------------------\n",
-      "|         211 |          26 |   1.697E+05 |   1.699E+05 |   1.639E+02 |   9.650E-04 |       40.27 |     3559.73 |\n",
+      "|         461 |           0 |   1.697E+05 |   1.699E+05 |   1.559E+02 |   9.176E-04 |       34.81 |     3565.19 |\n",
       "-----------------------------------------------------------------------------------------------------------------\n",
       " \n",
-      "Relative Tolerance Achieved\n",
-      "Optimal Solution Found at Node 1\n",
-      "Lower Bound: 169707.7780715559\n",
-      "Upper Bound: 169871.6980228018\n",
+      "Empty Stack: Exhaustive Search Finished\n",
+      "Optimal Solution Found at Node 71\n",
+      "Lower Bound: 169714.1201098805\n",
+      "Upper Bound: 169869.9993163163\n",
       "Solution:\n",
-      "   x[1] = 95.02153395384556\n",
-      "   x[2] = 0.2630982574134902\n",
-      "   x[3] = 0.01385685924905129\n",
-      "   x[4] = 0.050032869643406176\n",
-      "   x[5] = 26.31670002164094\n",
-      "   x[6] = 8.45937927400592\n",
+      "   x[1] = 68.70483394725427\n",
+      "   x[2] = 0.013856859244147913\n",
+      "   x[3] = 0.26309825737604337\n",
+      "   x[4] = 0.9499671303666462\n",
+      "   x[5] = 26.316700021840948\n",
+      "   x[6] = 8.459379273521812\n",
       " \n",
       "Termination Status: OPTIMAL\n",
       "Primal Status: FEASIBLE_POINT\n",
-      "Solve Time: 40.272\n",
-      "f^* = 169871.69802\n",
-      "x* = [95.022, 0.263, 0.014, 0.05, 26.317, 8.459]\n"
+      "Solve Time: 34.81\n",
+      "f^* = 169869.99932\n",
+      "x* = [68.705, 0.014, 0.263, 0.95, 26.317, 8.459]\n"
      ]
     }
    ],
    "source": [
-    "JuMP.optimize!(model)\n",
-    "println(\"Termination Status: $(JuMP.termination_status(model))\")\n",
-    "println(\"Primal Status: $(JuMP.primal_status(model))\")\n",
-    "println(\"Solve Time: $(round.(JuMP.solve_time(model), digits=5))\")\n",
-    "println(\"f^* = $(round(JuMP.objective_value(model), digits=5))\")\n",
-    "println(\"x* = $(round.(JuMP.value.(x), digits=3))\")"
+    "optimize!(model)\n",
+    "println(\"Termination Status: $(termination_status(model))\")\n",
+    "println(\"Primal Status: $(primal_status(model))\")\n",
+    "println(\"Solve Time: $(round.(solve_time(model), digits=5))\")\n",
+    "println(\"f^* = $(round(objective_value(model), digits=5))\")\n",
+    "println(\"x* = $(round.(value.(x), digits=3))\")"
    ]
   },
   {
@@ -464,26 +455,26 @@
      "data": {
       "text/plain": [
        "Dict{Any, Any} with 44 entries:\n",
-       "  cstr₊in₊F(t)        => 95.0215\n",
-       "  sep2₊outV₊F(t)      => 25.0\n",
-       "  sep2₊in₊y_C(t)      => 0.0500329\n",
-       "  sep1₊outV₊y_C(t)    => 0.0\n",
-       "  influent₊out₊y_B(t) => 0.0\n",
-       "  mixer₊in2₊y_A(t)    => 1.0\n",
-       "  cstr₊out₊y_A(t)     => 0.723045\n",
-       "  mixer₊in2₊F(t)      => 68.7048\n",
-       "  cstr₊in₊y_B(t)      => 0.0\n",
-       "  sep2₊outV₊y_A(t)    => 0.0\n",
-       "  mixer₊out₊y_C(t)    => 0.0\n",
-       "  sep1₊outV₊y_A(t)    => 1.0\n",
-       "  mixer₊in2₊y_B(t)    => -0.0\n",
-       "  mixer₊out₊y_B(t)    => 0.0\n",
-       "  mixer₊in1₊y_A(t)    => 1.0\n",
-       "  cstr₊in₊y_C(t)      => 0.0\n",
-       "  sep2₊in₊y_B(t)      => 0.949967\n",
-       "  sep2₊in₊y_A(t)      => 0.0\n",
-       "  mixer₊out₊F(t)      => 95.0215\n",
-       "  ⋮                   => ⋮"
+       "  sep1₊outV₊F(t)   => \u001b[34m68.7048\u001b[39m\u001b[0m\n",
+       "  sep1₊outV₊y_B(t) => \u001b[34m0\u001b[39m\u001b[0m\n",
+       "  mixer₊in2₊y_B(t) => \u001b[34m0\u001b[39m\u001b[0m\n",
+       "  sep1₊outL₊F(t)   => \u001b[34m26.3167\u001b[39m\u001b[0m\n",
+       "  sep2₊in₊y_B(t)   => \u001b[34m0.949967\u001b[39m\u001b[0m\n",
+       "  sep2₊outV₊y_A(t) => \u001b[34m0\u001b[39m\u001b[0m\n",
+       "  sep1₊in₊y_B(t)   => \u001b[34m0.263098\u001b[39m\u001b[0m\n",
+       "  sep2₊outL₊F(t)   => \u001b[34m1.3167\u001b[39m\u001b[0m\n",
+       "  sep2₊outV₊y_B(t) => \u001b[34m1.0\u001b[39m\u001b[0m\n",
+       "  mixer₊in1₊y_A(t) => \u001b[34m1.0\u001b[39m\u001b[0m\n",
+       "  sep1₊in₊y_C(t)   => \u001b[34m0.0138569\u001b[39m\u001b[0m\n",
+       "  sep1₊outV₊y_C(t) => \u001b[34m0\u001b[39m\u001b[0m\n",
+       "  sep2₊outV₊y_C(t) => \u001b[34m0\u001b[39m\u001b[0m\n",
+       "  mixer₊out₊y_B(t) => \u001b[34m0.0\u001b[39m\u001b[0m\n",
+       "  sep1₊outL₊y_A(t) => \u001b[34m0\u001b[39m\u001b[0m\n",
+       "  sep2₊in₊F(t)     => \u001b[34m26.3167\u001b[39m\u001b[0m\n",
+       "  cstr₊in₊y_C(t)   => \u001b[34m0.0\u001b[39m\u001b[0m\n",
+       "  sep2₊in₊y_A(t)   => \u001b[34m0\u001b[39m\u001b[0m\n",
+       "  mixer₊out₊F(t)   => \u001b[34m95.0215\u001b[39m\u001b[0m\n",
+       "  ⋮                => ⋮"
       ]
      },
      "metadata": {},
@@ -504,7 +495,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 11,
+   "execution_count": null,
    "id": "9c6d12f3",
    "metadata": {},
    "outputs": [
@@ -515,96 +506,19 @@
       "-----------------------------------------------------------------------------------------------------------------\n",
       "| Iteration # |    Nodes    | Lower Bound | Upper Bound |     Gap     |    Ratio    |    Timer    |  Time Left  |\n",
       "-----------------------------------------------------------------------------------------------------------------\n",
-      "|        1000 |         131 |   1.262E+05 |   1.699E+05 |   4.370E+04 |   2.572E-01 |       28.14 |     3571.86 |\n",
-      "|        2000 |         323 |   1.438E+05 |   1.699E+05 |   2.612E+04 |   1.537E-01 |       37.29 |     3562.71 |\n",
-      "|        3000 |         445 |   1.497E+05 |   1.699E+05 |   2.017E+04 |   1.187E-01 |       45.34 |     3554.66 |\n",
-      "|        4000 |         505 |   1.524E+05 |   1.699E+05 |   1.746E+04 |   1.028E-01 |       53.10 |     3546.90 |\n",
-      "|        5000 |         555 |   1.542E+05 |   1.699E+05 |   1.565E+04 |   9.211E-02 |       60.99 |     3539.01 |\n",
-      "|        6000 |         639 |   1.555E+05 |   1.699E+05 |   1.437E+04 |   8.458E-02 |       68.59 |     3531.41 |\n",
-      "|        7000 |         679 |   1.565E+05 |   1.699E+05 |   1.336E+04 |   7.863E-02 |       75.86 |     3524.14 |\n",
-      "|        8000 |         749 |   1.573E+05 |   1.699E+05 |   1.257E+04 |   7.401E-02 |       83.20 |     3516.80 |\n",
-      "|        9000 |         791 |   1.582E+05 |   1.699E+05 |   1.163E+04 |   6.847E-02 |       90.36 |     3509.64 |\n",
-      "|       10000 |         833 |   1.591E+05 |   1.699E+05 |   1.072E+04 |   6.313E-02 |       97.41 |     3502.59 |\n",
-      "|       11000 |         821 |   1.601E+05 |   1.699E+05 |   9.815E+03 |   5.778E-02 |      104.70 |     3495.30 |\n",
-      "|       12000 |         733 |   1.617E+05 |   1.699E+05 |   8.218E+03 |   4.838E-02 |      112.22 |     3487.78 |\n",
-      "|       13000 |         743 |   1.631E+05 |   1.699E+05 |   6.749E+03 |   3.973E-02 |      119.31 |     3480.69 |\n",
-      "|       14000 |         825 |   1.642E+05 |   1.699E+05 |   5.688E+03 |   3.348E-02 |      126.75 |     3473.25 |\n",
-      "|       15000 |         853 |   1.649E+05 |   1.699E+05 |   4.946E+03 |   2.912E-02 |      134.36 |     3465.64 |\n",
-      "|       16000 |         899 |   1.655E+05 |   1.699E+05 |   4.363E+03 |   2.569E-02 |      141.56 |     3458.44 |\n",
-      "|       17000 |         845 |   1.660E+05 |   1.699E+05 |   3.867E+03 |   2.276E-02 |      148.17 |     3451.83 |\n",
-      "|       18000 |         763 |   1.665E+05 |   1.699E+05 |   3.337E+03 |   1.964E-02 |      154.60 |     3445.40 |\n",
-      "|       19000 |         691 |   1.670E+05 |   1.699E+05 |   2.838E+03 |   1.671E-02 |      161.39 |     3438.61 |\n",
-      "|       20000 |         623 |   1.675E+05 |   1.699E+05 |   2.388E+03 |   1.406E-02 |      168.09 |     3431.91 |\n",
-      "|       21000 |         527 |   1.681E+05 |   1.699E+05 |   1.756E+03 |   1.034E-02 |      174.54 |     3425.46 |\n",
-      "|       22000 |         505 |   1.686E+05 |   1.699E+05 |   1.291E+03 |   7.601E-03 |      181.83 |     3418.17 |\n",
-      "|       23000 |         477 |   1.689E+05 |   1.699E+05 |   9.269E+02 |   5.456E-03 |      188.51 |     3411.49 |\n",
-      "|       24000 |         467 |   1.694E+05 |   1.699E+05 |   5.148E+02 |   3.031E-03 |      195.35 |     3404.65 |\n",
-      "|       25000 |         441 |   1.695E+05 |   1.699E+05 |   3.263E+02 |   1.921E-03 |      202.97 |     3397.03 |\n",
-      "|       26000 |         499 |   1.696E+05 |   1.699E+05 |   2.226E+02 |   1.311E-03 |      209.12 |     3390.88 |\n",
-      "|       26491 |         568 |   1.697E+05 |   1.699E+05 |   1.698E+02 |   9.995E-04 |      212.11 |     3387.89 |\n",
-      "-----------------------------------------------------------------------------------------------------------------\n",
-      " \n",
-      "Relative Tolerance Achieved\n",
-      "Optimal Solution Found at Node 1\n",
-      "Lower Bound: 169701.91518583894\n",
-      "Upper Bound: 169871.6980228018\n",
-      "Solution:\n",
-      "    x[1] = 26.31670002196759\n",
-      "    x[2] = 1.0\n",
-      "    x[3] = -2.990747056021119e-45\n",
-      "    x[4] = -3.905569389055896e-44\n",
-      "    x[5] = 26.31670002196759\n",
-      "    x[6] = 1.0\n",
-      "    x[7] = 1.1976322549267912e-44\n",
-      "    x[8] = -2.912116506157164e-44\n",
-      "    x[9] = 68.7048339160128\n",
-      "   x[10] = 1.0\n",
-      "   x[11] = -7.364466832721216e-45\n",
-      "   x[12] = -4.1014804181443646e-43\n",
-      "   x[13] = 95.0215339379804\n",
-      "   x[14] = 1.0\n",
-      "   x[15] = -5.114739394785582e-44\n",
-      "   x[16] = -4.023413752530357e-43\n",
-      "   x[17] = 95.0215339379804\n",
-      "   x[18] = 1.0\n",
-      "   x[19] = -7.496946784137771e-44\n",
-      "   x[20] = -3.914877584707458e-43\n",
-      "   x[21] = 95.0215339379804\n",
-      "   x[22] = 0.7230448832877795\n",
-      "   x[23] = 0.2630982574574181\n",
-      "   x[24] = 0.01385685925480256\n",
-      "   x[25] = 95.0215339379804\n",
-      "   x[26] = 0.7230448832877795\n",
-      "   x[27] = 0.2630982574574181\n",
-      "   x[28] = 0.01385685925480256\n",
-      "   x[29] = 68.7048339160128\n",
-      "   x[30] = 1.0\n",
-      "   x[31] = 0.0\n",
-      "   x[32] = 1.4743240170669853e-45\n",
-      "   x[33] = 26.316700021967595\n",
-      "   x[34] = -3.4035930753771503e-53\n",
-      "   x[35] = 0.9499671303448025\n",
-      "   x[36] = 0.05003286965519755\n",
-      "   x[37] = 26.316700021967595\n",
-      "   x[38] = -4.176194859519056e-53\n",
-      "   x[39] = 0.9499671303448025\n",
-      "   x[40] = 0.05003286965519755\n",
-      "   x[41] = 25.000000000013557\n",
-      "   x[42] = -3.330793643180377e-60\n",
-      "   x[43] = 1.0\n",
-      "   x[44] = 0.0\n",
-      "   x[45] = 1.3167000219540375\n",
-      "   x[46] = -6.372367644529809e-58\n",
-      "   x[47] = 1.6313261169996311e-55\n",
-      "   x[48] = 1.0\n",
-      "   x[49] = 26.31670002196759\n",
-      "   x[50] = 8.459379274754202\n",
-      " \n",
-      "Termination Status: OPTIMAL\n",
-      "Primal Status: FEASIBLE_POINT\n",
-      "Solve Time: 212.113\n",
-      "f^* = 169871.69802\n",
-      "x* = [26.317, 1.0, -0.0, -0.0, 26.317, 1.0, 0.0, -0.0, 68.705, 1.0, -0.0, -0.0, 95.022, 1.0, -0.0, -0.0, 95.022, 1.0, -0.0, -0.0, 95.022, 0.723, 0.263, 0.014, 95.022, 0.723, 0.263, 0.014, 68.705, 1.0, 0.0, 0.0, 26.317, -0.0, 0.95, 0.05, 26.317, -0.0, 0.95, 0.05, 25.0, -0.0, 1.0, 0.0, 1.317, -0.0, 0.0, 1.0, 26.317, 8.459]\n"
+      "|        1000 |         159 |   1.031E+05 |   1.699E+05 |   6.677E+04 |   3.930E-01 |       24.43 |     3575.57 |\n",
+      "|        2000 |         251 |   1.244E+05 |   1.699E+05 |   4.552E+04 |   2.680E-01 |       25.09 |     3574.91 |\n",
+      "|        3000 |         265 |   1.358E+05 |   1.699E+05 |   3.411E+04 |   2.008E-01 |       25.78 |     3574.22 |\n",
+      "|        4000 |         349 |   1.419E+05 |   1.699E+05 |   2.800E+04 |   1.648E-01 |       26.51 |     3573.49 |\n",
+      "|        5000 |         341 |   1.473E+05 |   1.699E+05 |   2.262E+04 |   1.332E-01 |       27.28 |     3572.72 |\n",
+      "|        6000 |         323 |   1.512E+05 |   1.699E+05 |   1.864E+04 |   1.098E-01 |       27.98 |     3572.02 |\n",
+      "|        7000 |         301 |   1.556E+05 |   1.699E+05 |   1.427E+04 |   8.400E-02 |       28.67 |     3571.33 |\n",
+      "|        8000 |         407 |   1.599E+05 |   1.699E+05 |   9.928E+03 |   5.844E-02 |       29.41 |     3570.59 |\n",
+      "|        9000 |         441 |   1.633E+05 |   1.699E+05 |   6.576E+03 |   3.871E-02 |       30.18 |     3569.82 |\n",
+      "|       10000 |         415 |   1.659E+05 |   1.699E+05 |   3.980E+03 |   2.343E-02 |       30.92 |     3569.08 |\n",
+      "|       11000 |         445 |   1.672E+05 |   1.699E+05 |   2.649E+03 |   1.559E-02 |       31.61 |     3568.39 |\n",
+      "|       12000 |         503 |   1.680E+05 |   1.699E+05 |   1.858E+03 |   1.094E-02 |       32.26 |     3567.74 |\n",
+      "|       13000 |         497 |   1.685E+05 |   1.699E+05 |   1.370E+03 |   8.066E-03 |       32.92 |     3567.08 |\n"
      ]
     }
    ],
@@ -617,12 +531,12 @@
     "\n",
     "register_nlsystem(full_model, full_system, obj, [g1, g2])\n",
     "\n",
-    "JuMP.optimize!(full_model)\n",
-    "println(\"Termination Status: $(JuMP.termination_status(full_model))\")\n",
-    "println(\"Primal Status: $(JuMP.primal_status(full_model))\")\n",
-    "println(\"Solve Time: $(round.(JuMP.solve_time(full_model), digits=5))\")\n",
-    "println(\"f^* = $(round(JuMP.objective_value(full_model), digits=5))\")\n",
-    "println(\"x* = $(round.(JuMP.value.(x), digits=3))\")"
+    "optimize!(full_model)\n",
+    "println(\"Termination Status: $(termination_status(full_model))\")\n",
+    "println(\"Primal Status: $(primal_status(full_model))\")\n",
+    "println(\"Solve Time: $(round.(solve_time(full_model), digits=5))\")\n",
+    "println(\"f^* = $(round(objective_value(full_model), digits=5))\")\n",
+    "println(\"x* = $(round.(value.(x), digits=3))\")"
    ]
   },
   {
@@ -641,7 +555,9 @@
     "### References\n",
     "\n",
     "1. Kokossis, A.C. and Floudas, C.A. Synthesis of isothermal reactor-separator-recycle systems. *Chemical Engineering Science.* 46, 1361-1383 (1991). DOI: [10.1016/0009-2509(91)85063-4](https://doi.org/10.1016/0009-2509(91)85063-4)\n",
-    "2. Wilhelm, M. E. and Stuber, M.D. EAGO.jl: easy advanced global optimization in Julia. *Optimization Methods and Software.* 37(2), 425-450 (2022). DOI: [10.1080/10556788.2020.1786566](https://doi.org/10.1080/10556788.2020.1786566)"
+    "2. Ma, Y., Gowda, S., Anantharaman, R., Laughman, C., Shah, V., and Rackauckas, C. ModelingToolkit: A Composable Graph Transformation System For Equation-Based Modeling. (2022). DOI: [10.48550/arXiv.2103.05244](https://doi.org/10.48550/arXiv.2103.05244)\n",
+    "3. Lubin, M., Dowson, O., Dias Garcia, J., Huchette, J., Legat, B., and Vielma, J.P. JuMP 1.0: recent improvements to a modeling language for mathematical optimization. *Mathematical Programming Computation.* 15, 581-589 (2023). DOI: [10.1007/s12532-023-00239-3](https://doi.org/10.1007/s12532-023-00239-3)\n",
+    "4. Wilhelm, M. E. and Stuber, M.D. EAGO.jl: easy advanced global optimization in Julia. *Optimization Methods & Software.* 37(2), 425-450 (2022). DOI: [10.1080/10556788.2020.1786566](https://doi.org/10.1080/10556788.2020.1786566)"
    ]
   }
  ],
diff --git a/examples/algebraic_model.jl b/examples/algebraic_model.jl
index c09ba1c..3c8beff 100644
--- a/examples/algebraic_model.jl
+++ b/examples/algebraic_model.jl
@@ -1,4 +1,7 @@
-using ModelingToolkit, JuMP, EAGO, EOptInterface
+using EAGO
+using EOptInterface
+using JuMP
+using ModelingToolkit
 using ModelingToolkit: t_nounits as t, D_nounits as D
 
 # Create algebraic ModelingToolkit system
@@ -163,14 +166,14 @@ xU = [100.0, 1.0, 1.0, 1.0, 100.0, 10.0]
 register_nlsystem(model, system, obj, [g1, g2])
 
 # Optimize model
-JuMP.optimize!(model)
+optimize!(model)
 
 # Display results
-println("Termination Status: $(JuMP.termination_status(model))")
-println("Primal Status: $(JuMP.primal_status(model))")
-println("Solve Time: $(round.(JuMP.solve_time(model), digits=5))")
-println("f^* = $(round(JuMP.objective_value(model), digits=5))")
-println("x* = $(round.(JuMP.value.(x), digits=3))")
+println("Termination Status: $(termination_status(model))")
+println("Primal Status: $(primal_status(model))")
+println("Solve Time: $(round.(solve_time(model), digits=5))")
+println("f^* = $(round(objective_value(model), digits=5))")
+println("x* = $(round.(value.(x), digits=3))")
 
 # Retrieve full-space solution
 full_solution(model, system)
@@ -191,9 +194,9 @@ xU = vcat(repeat([100.0, 1.0, 1.0, 1.0], 12), 100.0, 10.0)
 register_nlsystem(full_model, full_system, obj, [g1, g2])
 
 # Optimize model and retrieve results
-JuMP.optimize!(full_model)
-println("Termination Status: $(JuMP.termination_status(full_model))")
-println("Primal Status: $(JuMP.primal_status(full_model))")
-println("Solve Time: $(round.(JuMP.solve_time(full_model), digits=5))")
-println("f^* = $(round(JuMP.objective_value(full_model), digits=5))")
-println("x* = $(round.(JuMP.value.(x), digits=3))")
+optimize!(full_model)
+println("Termination Status: $(termination_status(full_model))")
+println("Primal Status: $(primal_status(full_model))")
+println("Solve Time: $(round.(solve_time(full_model), digits=5))")
+println("f^* = $(round(objective_value(full_model), digits=5))")
+println("x* = $(round.(value.(x), digits=3))")
diff --git a/examples/kinetic_intensity_data.csv b/examples/kinetic_intensity_data.csv
new file mode 100644
index 0000000..7970fab
--- /dev/null
+++ b/examples/kinetic_intensity_data.csv
@@ -0,0 +1,202 @@
+time,intensity
+0,0
+0.01,66.0952
+0.02,104.762
+0.03,110.333
+0.04,114.905
+0.05,122.238
+0.06,125.429
+0.07,125.429
+0.08,123.476
+0.09,121.286
+0.1,118.857
+0.11,117.667
+0.12,116.143
+0.13,113.857
+0.14,111.571
+0.15,108.81
+0.16,105.952
+0.17,104.048
+0.18,102.048
+0.19,100.143
+0.2,98.5238
+0.21,96.2381
+0.22,94.381
+0.23,91.6667
+0.24,89.5714
+0.25,87.1429
+0.26,84.8571
+0.27,83.4286
+0.28,81.1905
+0.29,78.9048
+0.3,77.0476
+0.31,75.4762
+0.32,73.4762
+0.33,71.8095
+0.34,70.6667
+0.35,68.381
+0.36,67.3333
+0.37,65.0952
+0.38,63.7143
+0.39,62.0476
+0.4,60.8571
+0.41,59.619
+0.42,58.2857
+0.43,57.4762
+0.44,56.4762
+0.45,55.8095
+0.46,54.5238
+0.47,53
+0.48,51.8571
+0.49,50.4286
+0.5,49.381
+0.51,47.9524
+0.52,47.3714
+0.53,46.8952
+0.54,46.4857
+0.55,45.9048
+0.56,45.0762
+0.57,44.3238
+0.58,43.4143
+0.59,43.5429
+0.6,42.3619
+0.61,41.8381
+0.62,40.2381
+0.63,39.1286
+0.64,38.7857
+0.65,37.081
+0.66,36.9524
+0.67,36.581
+0.68,36.281
+0.69,35.3476
+0.7,34.8905
+0.71,34.1667
+0.72,33.6714
+0.73,32.9667
+0.74,31.8429
+0.75,31.5429
+0.76,31.1476
+0.77,30.9905
+0.78,29.9571
+0.79,29.1333
+0.8,28.7857
+0.81,28.4429
+0.82,28.3476
+0.83,27.5429
+0.84,27.4333
+0.85,27.6048
+0.86,27.1762
+0.87,27.2
+0.88,26.4333
+0.89,25.7619
+0.9,24.8095
+0.91,24.7429
+0.92,24.2857
+0.93,24.1714
+0.94,23.5667
+0.95,23.5476
+0.96,23.3952
+0.97,22.919
+0.98,22.3095
+0.99,21.8048
+1,21.2857
+1.01,21.2048
+1.02,20.8429
+1.03,20.4429
+1.04,20.0048
+1.05,19.9381
+1.06,19.5
+1.07,19.8667
+1.08,18.9333
+1.09,19.1381
+1.1,18.9619
+1.11,18.5476
+1.12,17.9048
+1.13,17.7571
+1.14,18.5333
+1.15,18.3762
+1.16,18.3571
+1.17,18.3286
+1.18,18.2762
+1.19,18.3952
+1.2,17.5952
+1.21,18.1524
+1.22,18.1952
+1.23,17.8476
+1.24,17.9095
+1.25,17.5048
+1.26,17.5
+1.27,15.9619
+1.28,16.2095
+1.29,16.181
+1.3,15.6952
+1.31,15.7095
+1.32,15.4619
+1.33,15.9476
+1.34,16
+1.35,16.1952
+1.36,16.1143
+1.37,15.7429
+1.38,15.5762
+1.39,15.7048
+1.4,15.8095
+1.41,15.6667
+1.42,14.9048
+1.43,14.5857
+1.44,14.7524
+1.45,14.7571
+1.46,14.9762
+1.47,14.5333
+1.48,14.5524
+1.49,14.0143
+1.5,13.6286
+1.51,13.4429
+1.52,13.4667
+1.53,13.319
+1.54,12.9333
+1.55,13.1238
+1.56,12.7476
+1.57,12.9333
+1.58,13.0714
+1.59,13.0714
+1.6,12.7619
+1.61,12.4238
+1.62,12.5143
+1.63,12.9143
+1.64,12.5714
+1.65,13.3667
+1.66,13.2286
+1.67,13.7905
+1.68,13.7571
+1.69,13.5905
+1.7,12.9667
+1.71,12.981
+1.72,12.8857
+1.73,12.919
+1.74,13.0143
+1.75,13.0095
+1.76,12.3857
+1.77,12.5571
+1.78,12.3429
+1.79,12.7571
+1.8,12.681
+1.81,12.5429
+1.82,12.1857
+1.83,12.7905
+1.84,12.5571
+1.85,12.8429
+1.86,12.5476
+1.87,12.5714
+1.88,12.3762
+1.89,11.9952
+1.9,11.4571
+1.91,11.3
+1.92,11.1524
+1.93,11.681
+1.94,11.619
+1.95,11.9048
+1.96,12
+1.97,12.0762
+1.98,11.9143
+1.99,11.7619
+2,11.5333
diff --git a/examples/kinetic_intensity_data.jl b/examples/kinetic_intensity_data.jl
deleted file mode 100644
index 353015b..0000000
--- a/examples/kinetic_intensity_data.jl
+++ /dev/null
@@ -1,202 +0,0 @@
-data = [
-    66.0952
-    104.762
-    110.333
-    114.905
-    122.238
-    125.429
-    125.429
-    123.476
-    121.286
-    118.857
-    117.667
-    116.143
-    113.857
-    111.571
-    108.81
-    105.952
-    104.048
-    102.048
-    100.143
-    98.5238
-    96.2381
-    94.381
-    91.6667
-    89.5714
-    87.1429
-    84.8571
-    83.4286
-    81.1905
-    78.9048
-    77.0476
-    75.4762
-    73.4762
-    71.8095
-    70.6667
-    68.381
-    67.3333
-    65.0952
-    63.7143
-    62.0476
-    60.8571
-    59.619
-    58.2857
-    57.4762
-    56.4762
-    55.8095
-    54.5238
-    53.0
-    51.8571
-    50.4286
-    49.381
-    47.9524
-    47.3714
-    46.8952
-    46.4857
-    45.9048
-    45.0762
-    44.3238
-    43.4143
-    43.5429
-    42.3619
-    41.8381
-    40.2381
-    39.1286
-    38.7857
-    37.081
-    36.9524
-    36.581
-    36.281
-    35.3476
-    34.8905
-    34.1667
-    33.6714
-    32.9667
-    31.8429
-    31.5429
-    31.1476
-    30.9905
-    29.9571
-    29.1333
-    28.7857
-    28.4429
-    28.3476
-    27.5429
-    27.4333
-    27.6048
-    27.1762
-    27.2
-    26.4333
-    25.7619
-    24.8095
-    24.7429
-    24.2857
-    24.1714
-    23.5667
-    23.5476
-    23.3952
-    22.919
-    22.3095
-    21.8048
-    21.2857
-    21.2048
-    20.8429
-    20.4429
-    20.0048
-    19.9381
-    19.5
-    19.8667
-    18.9333
-    19.1381
-    18.9619
-    18.5476
-    17.9048
-    17.7571
-    18.5333
-    18.3762
-    18.3571
-    18.3286
-    18.2762
-    18.3952
-    17.5952
-    18.1524
-    18.1952
-    17.8476
-    17.9095
-    17.5048
-    17.5
-    15.9619
-    16.2095
-    16.181
-    15.6952
-    15.7095
-    15.4619
-    15.9476
-    16.0
-    16.1952
-    16.1143
-    15.7429
-    15.5762
-    15.7048
-    15.8095
-    15.6667
-    14.9048
-    14.5857
-    14.7524
-    14.7571
-    14.9762
-    14.5333
-    14.5524
-    14.0143
-    13.6286
-    13.4429
-    13.4667
-    13.319
-    12.9333
-    13.1238
-    12.7476
-    12.9333
-    13.0714
-    13.0714
-    12.7619
-    12.4238
-    12.5143
-    12.9143
-    12.5714
-    13.3667
-    13.2286
-    13.7905
-    13.7571
-    13.5905
-    12.9667
-    12.981
-    12.8857
-    12.919
-    13.0143
-    13.0095
-    12.3857
-    12.5571
-    12.3429
-    12.7571
-    12.681
-    12.5429
-    12.1857
-    12.7905
-    12.5571
-    12.8429
-    12.5476
-    12.5714
-    12.3762
-    11.9952
-    11.4571
-    11.3
-    11.1524
-    11.681
-    11.619
-    11.9048
-    12.0
-    12.0762
-    11.9143
-    11.7619
-    11.5333
-]
\ No newline at end of file
diff --git a/examples/ode_model.ipynb b/examples/ode_model.ipynb
index cf0a133..47da60c 100644
--- a/examples/ode_model.ipynb
+++ b/examples/ode_model.ipynb
@@ -7,7 +7,7 @@
    "source": [
     "# Kinetic Parameter Estimation\n",
     "\n",
-    "A kinetic parameter estimation problem originally described by [[1](#references)] is used to demonstrate the use of EOptInterface to formulate a `JuMP.Model` from an ODE `ModelingToolkit.System`. The reacting system of interest is represented by the following initial value problem:\n",
+    "A kinetic parameter estimation problem originally described by [[1](#references)] is used to demonstrate the use of EOptInterface to formulate a `JuMP.Model` [[2](#references)] from an ODE `ModelingToolkit.System` [[3](#references)]. The reacting system of interest is represented by the following initial value problem:\n",
     "$$\n",
     "\\begin{align*}\n",
     "    &\\frac{dx_A}{dt} = k_1 x_Z x_Y - c_{\\text{O}_{2}} (k_{2f} + k_{3f}) x_A + \\frac{k_{2f}}{K_2} x_D + \\frac{k_{3f}}{K_3} x_B - k_5 x_A^2 \\\\\n",
@@ -23,12 +23,17 @@
   },
   {
    "cell_type": "code",
-   "execution_count": null,
+   "execution_count": 1,
    "id": "341eed5e",
    "metadata": {},
    "outputs": [],
    "source": [
-    "using ModelingToolkit, JuMP, Ipopt, EOptInterface\n",
+    "using CSV\n",
+    "using DataFrames\n",
+    "using EOptInterface\n",
+    "using Ipopt\n",
+    "using JuMP\n",
+    "using ModelingToolkit\n",
     "using ModelingToolkit: t_nounits as t, D_nounits as D"
    ]
   },
@@ -42,7 +47,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 11,
+   "execution_count": 2,
    "id": "28356dee",
    "metadata": {},
    "outputs": [],
@@ -91,28 +96,26 @@
    "id": "a16811dc",
    "metadata": {},
    "source": [
-    "The objective is to determine the values of uncertain kinetic parameters that best fit this model against experimental intensity versus time data available in `kinetic_intensity_data.jl`, where the relationship between intensity and concentration is known to be $I^{\\text{calc}} = x_A + \\frac{2}{21} x_B + \\frac{2}{21} x_D$ [[2](#references)]. The experimental data set $I^{\\text{exp}}$ contains 200 measurements between $(0.01, 2.00)$ spaced linearly by $\\Delta t = 0.01$. Based on this information, we define the integration time span to be $(0, 2)$ to include the initial condition, and use the same time step as between measurements, resulting in $N = 201$ discrete time points. The objective function is then defined as sum of squared error between the model-predicted and measured intensity:\n",
+    "The objective is to determine the values of uncertain kinetic parameters that best fit this model against experimental intensity versus time data available in `kinetic_intensity_data.csv`, where the relationship between intensity and concentration is known to be $I^{\\text{calc}} = x_A + \\frac{2}{21} x_B + \\frac{2}{21} x_D$ [[4](#references)]. The experimental data set $I^{\\text{exp}}$ contains 201 measurements between $(0, 2)$ spaced linearly by $\\Delta t = 0.01$. Based on this information, we define the integration time span to be $(0, 2)$ to include the initial condition, and use the same time step as between measurements, resulting in $N = 201$ discrete time points. The objective function is then defined as sum of squared error between the model-predicted and measured intensity:\n",
     "$$\n",
     "\\begin{align}\n",
-    "    f(\\mathbf{x}) = \\sum_{i=2}^{N} \\left(I^{\\text{calc}}(\\mathbf{x}_{i}) - I^{\\text{exp}}_{i-1} \\right)^2,\n",
+    "    f(\\mathbf{x}) = \\sum_{i=2}^{N} \\left(I^{\\text{calc}}(\\mathbf{x}_{i}) - I^{\\text{exp}}_{i} \\right)^2.\n",
     "\\end{align}\n",
-    "$$\n",
-    "where the offset in indexing is due to the initial condition being stored at $i = 1$, i.e., $\\mathbf x_1 = (x_{Z,0}, x_{Y,0}, x_{D,0}, x_{B,0}, x_{A,0})$ and the experimental\n",
-    "data vector only containing measured data after the experiment begins (i.e., after the initial condition)."
+    "$$"
    ]
   },
   {
    "cell_type": "code",
-   "execution_count": 12,
+   "execution_count": 3,
    "id": "413b1871",
    "metadata": {},
    "outputs": [],
    "source": [
-    "tspan = (0.0, 2.0)\n",
-    "tstep = 0.01\n",
-    "N = Int(floor((tspan[2] - tspan[1])/tstep)) + 1\n",
-    "include(\"kinetic_intensity_data.jl\")\n",
-    "intensity(x_A, x_B, x_D) = x_A + 2/21*x_B + 2/21*x_D;"
+    "data = CSV.read(\"kinetic_intensity_data.csv\", DataFrame)\n",
+    "intensity(x_A, x_B, x_D) = x_A + 2/21*x_B + 2/21*x_D\n",
+    "tspan = (data.time[1], data.time[end])\n",
+    "tstep = data.time[2] - data.time[1]\n",
+    "N = length(data.time);"
    ]
   },
   {
@@ -120,12 +123,12 @@
    "id": "9f776dac",
    "metadata": {},
    "source": [
-    "We can then define the `JuMP.Model` using an appropriate solver of choice. Given that this is a nonlinear system, we have chosen to use `Ipopt` [[3](#references)]."
+    "We can then define the `JuMP.Model` using an appropriate solver of choice. Given that this is a nonlinear system, we have chosen to use `Ipopt` [[5](#references)]."
    ]
   },
   {
    "cell_type": "code",
-   "execution_count": 13,
+   "execution_count": 4,
    "id": "5ca7a46f",
    "metadata": {},
    "outputs": [],
@@ -143,14 +146,14 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 14,
+   "execution_count": 5,
    "id": "8f38e3b9",
    "metadata": {},
    "outputs": [
     {
      "data": {
       "text/plain": [
-       "8-element Vector{Any}:\n",
+       "8-element Vector{SymbolicUtils.BasicSymbolicImpl.var\"typeof(BasicSymbolicImpl)\"{SymReal}}:\n",
        " x_Z(t)\n",
        " x_Y(t)\n",
        " x_D(t)\n",
@@ -179,7 +182,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 15,
+   "execution_count": 6,
    "id": "7a7f6375",
    "metadata": {},
    "outputs": [],
@@ -201,7 +204,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 16,
+   "execution_count": 7,
    "id": "9b15d798",
    "metadata": {},
    "outputs": [],
@@ -219,12 +222,12 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 17,
+   "execution_count": 8,
    "id": "d9762428",
    "metadata": {},
    "outputs": [],
    "source": [
-    "@objective(model, Min, sum((intensity(z[5,i], z[4,i], z[3,i]) - data[i-1])^2 for i in 2:N));"
+    "@objective(model, Min, sum((intensity(z[5,i], z[4,i], z[3,i]) - data.intensity[i])^2 for i in 1:N));"
    ]
   },
   {
@@ -237,7 +240,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 18,
+   "execution_count": 9,
    "id": "a04a5a8e",
    "metadata": {},
    "outputs": [
@@ -245,6 +248,13 @@
      "name": "stdout",
      "output_type": "stream",
      "text": [
+      "\n",
+      "******************************************************************************\n",
+      "This program contains Ipopt, a library for large-scale nonlinear optimization.\n",
+      " Ipopt is released as open source code under the Eclipse Public License (EPL).\n",
+      "         For more information visit https://github.com/coin-or/Ipopt\n",
+      "******************************************************************************\n",
+      "\n",
       "This is Ipopt version 3.14.19, running with linear solver MUMPS 5.8.2.\n",
       "\n",
       "Number of nonzeros in equality constraint Jacobian...:     6979\n",
@@ -295,17 +305,17 @@
       "  28  9.6227604e+03 1.76e-04 7.62e-04  -3.8 3.22e-01    -  1.00e+00 1.00e+00h  1\n",
       "  29  9.6227628e+03 3.75e-07 4.30e-07  -5.7 1.11e-02    -  1.00e+00 1.00e+00h  1\n",
       "iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls\n",
-      "  30  9.6227629e+03 1.98e-13 2.62e-12  -8.6 2.62e-05    -  1.00e+00 1.00e+00h  1\n",
+      "  30  9.6227629e+03 1.98e-13 7.96e-13  -8.6 2.62e-05    -  1.00e+00 1.00e+00h  1\n",
       "\n",
       "Number of Iterations....: 30\n",
       "\n",
       "                                   (scaled)                 (unscaled)\n",
-      "Objective...............:   3.8359401942857030e+03    9.6227628525812288e+03\n",
-      "Dual infeasibility......:   2.6152565902065386e-12    6.5605803770603188e-12\n",
-      "Constraint violation....:   1.9795276529066541e-13    1.9795276529066541e-13\n",
+      "Objective...............:   3.8359401942857071e+03    9.6227628525812397e+03\n",
+      "Dual infeasibility......:   7.9626718666068208e-13    1.9974999391132541e-12\n",
+      "Constraint violation....:   1.9773072068574038e-13    1.9773072068574038e-13\n",
       "Variable bound violation:   0.0000000000000000e+00    0.0000000000000000e+00\n",
-      "Complementarity.........:   2.5066667739822905e-09    6.2881741358764948e-09\n",
-      "Overall NLP error.......:   2.5066667739822905e-09    6.2881741358764948e-09\n",
+      "Complementarity.........:   2.5066667739789243e-09    6.2881741358680501e-09\n",
+      "Overall NLP error.......:   2.5066667739789243e-09    6.2881741358680501e-09\n",
       "\n",
       "\n",
       "Number of objective function evaluations             = 32\n",
@@ -315,24 +325,24 @@
       "Number of equality constraint Jacobian evaluations   = 31\n",
       "Number of inequality constraint Jacobian evaluations = 0\n",
       "Number of Lagrangian Hessian evaluations             = 30\n",
-      "Total seconds in IPOPT                               = 0.157\n",
+      "Total seconds in IPOPT                               = 0.233\n",
       "\n",
       "EXIT: Optimal Solution Found.\n",
       "Termination Status: LOCALLY_SOLVED\n",
       "Primal Status: FEASIBLE_POINT\n",
-      "Solve Time: 0.159\n",
+      "Solve Time: 0.243\n",
       "f^* = 9622.76285\n",
       "p* = [828.065, 385.732, 14.567]\n"
      ]
     }
    ],
    "source": [
-    "JuMP.optimize!(model)\n",
-    "println(\"Termination Status: $(JuMP.termination_status(model))\")\n",
-    "println(\"Primal Status: $(JuMP.primal_status(model))\")\n",
-    "println(\"Solve Time: $(round.(JuMP.solve_time(model), digits=5))\")\n",
-    "println(\"f^* = $(round(JuMP.objective_value(model), digits=5))\")\n",
-    "println(\"p* = $(round.(JuMP.value.(p), digits=3))\")"
+    "optimize!(model)\n",
+    "println(\"Termination Status: $(termination_status(model))\")\n",
+    "println(\"Primal Status: $(primal_status(model))\")\n",
+    "println(\"Solve Time: $(round.(solve_time(model), digits=5))\")\n",
+    "println(\"f^* = $(round(objective_value(model), digits=5))\")\n",
+    "println(\"p* = $(round.(value.(p), digits=3))\")"
    ]
   },
   {
@@ -343,8 +353,10 @@
     "### References\n",
     "\n",
     "1. Taylor, J.W. Direct measurement and analysis of cyclohexadienyl oxidation. Ph.D. thesis, Massachusetts Institute of Technology. (2005). URL: http://hdl.handle.net/1721.1/33716\n",
-    "2. Singer, A.B. Global dynamic optimization. Ph.D. thesis, Massachusetts Institute of Technology. (2004). URL: http://hdl.handle.net/1721.1/28662\n",
-    "3. Wächter, A. and Biegler, L.T. On the implementation of an interior-point filter line-search algorithm for large-scale nonlinear programming. *Mathematical Programming.* 106, 25-57 (2006). DOI: [10.1007/s10107-004-0559-y](https://doi.org/10.1007/s10107-004-0559-y)"
+    "2. Lubin, M., Dowson, O., Dias Garcia, J., Huchette, J., Legat, B., and Vielma, J.P. JuMP 1.0: recent improvements to a modeling language for mathematical optimization. *Mathematical Programming Computation.* 15, 581-589 (2023). DOI: [10.1007/s12532-023-00239-3](https://doi.org/10.1007/s12532-023-00239-3)\n",
+    "3. Ma, Y., Gowda, S., Anantharaman, R., Laughman, C., Shah, V., and Rackauckas, C. ModelingToolkit: A Composable Graph Transformation System For Equation-Based Modeling. (2022). DOI: [10.48550/arXiv.2103.05244](https://doi.org/10.48550/arXiv.2103.05244)\n",
+    "4. Singer, A.B. Global dynamic optimization. Ph.D. thesis, Massachusetts Institute of Technology. (2004). URL: http://hdl.handle.net/1721.1/28662\n",
+    "5. Wächter, A. and Biegler, L.T. On the implementation of an interior-point filter line-search algorithm for large-scale nonlinear programming. *Mathematical Programming.* 106, 25-57 (2006). DOI: [10.1007/s10107-004-0559-y](https://doi.org/10.1007/s10107-004-0559-y)"
    ]
   }
  ],
diff --git a/examples/ode_model.jl b/examples/ode_model.jl
index 2c5d4a4..87c0c1d 100644
--- a/examples/ode_model.jl
+++ b/examples/ode_model.jl
@@ -1,4 +1,9 @@
-using ModelingToolkit, JuMP, Ipopt, EOptInterface
+using CSV
+using DataFrames
+using EOptInterface
+using Ipopt
+using JuMP
+using ModelingToolkit
 using ModelingToolkit: t_nounits as t, D_nounits as D
 
 # Create ModelingToolkit ODE system
@@ -41,15 +46,14 @@ end
 # Compile system
 @mtkcompile system = KineticParameterEstimation()
 
-# Define integration time span
-tspan = (0.0, 2.0)
-tstep = 0.01
-N = Int(floor((tspan[2] - tspan[1])/tstep)) + 1
-
 # Include experimental intensity data
-include("kinetic_intensity_data.jl")
+data = CSV.read("examples/kinetic_intensity_data.csv", DataFrame)
 # Define intensity function
 intensity(x_A, x_B, x_D) = x_A + 2/21*x_B + 2/21*x_D
+# Retrieve integration time data
+tspan = (data.time[1], data.time[end])
+tstep = data.time[2] - data.time[1]
+N = length(data.time)
 
 # Create JuMP model
 model = Model(Ipopt.Optimizer)
@@ -73,12 +77,12 @@ pU = [1200.0, 1200.0, 40.0]
 register_odesystem(model, system, tspan, tstep, "IE")
 
 # Define objective function
-@objective(model, Min, sum((intensity(z[5,i], z[4,i], z[3,i]) - data[i-1])^2 for i in 2:N))
+@objective(model, Min, sum((intensity(z[5,i], z[4,i], z[3,i]) - data.intensity[i])^2 for i in 1:N))
 
 # Optimize model and retrieve results
-JuMP.optimize!(model)
-println("Termination Status: $(JuMP.termination_status(model))")
-println("Primal Status: $(JuMP.primal_status(model))")
-println("Solve Time: $(round.(JuMP.solve_time(model), digits=5))")
-println("f^* = $(round(JuMP.objective_value(model), digits=5))")
-println("p* = $(round.(JuMP.value.(p), digits=3))")
+optimize!(model)
+println("Termination Status: $(termination_status(model))")
+println("Primal Status: $(primal_status(model))")
+println("Solve Time: $(round.(solve_time(model), digits=5))")
+println("f^* = $(round(objective_value(model), digits=5))")
+println("p* = $(round.(value.(p), digits=3))")
diff --git a/test/kinetic_intensity_data.csv b/test/kinetic_intensity_data.csv
new file mode 100644
index 0000000..7970fab
--- /dev/null
+++ b/test/kinetic_intensity_data.csv
@@ -0,0 +1,202 @@
+time,intensity
+0,0
+0.01,66.0952
+0.02,104.762
+0.03,110.333
+0.04,114.905
+0.05,122.238
+0.06,125.429
+0.07,125.429
+0.08,123.476
+0.09,121.286
+0.1,118.857
+0.11,117.667
+0.12,116.143
+0.13,113.857
+0.14,111.571
+0.15,108.81
+0.16,105.952
+0.17,104.048
+0.18,102.048
+0.19,100.143
+0.2,98.5238
+0.21,96.2381
+0.22,94.381
+0.23,91.6667
+0.24,89.5714
+0.25,87.1429
+0.26,84.8571
+0.27,83.4286
+0.28,81.1905
+0.29,78.9048
+0.3,77.0476
+0.31,75.4762
+0.32,73.4762
+0.33,71.8095
+0.34,70.6667
+0.35,68.381
+0.36,67.3333
+0.37,65.0952
+0.38,63.7143
+0.39,62.0476
+0.4,60.8571
+0.41,59.619
+0.42,58.2857
+0.43,57.4762
+0.44,56.4762
+0.45,55.8095
+0.46,54.5238
+0.47,53
+0.48,51.8571
+0.49,50.4286
+0.5,49.381
+0.51,47.9524
+0.52,47.3714
+0.53,46.8952
+0.54,46.4857
+0.55,45.9048
+0.56,45.0762
+0.57,44.3238
+0.58,43.4143
+0.59,43.5429
+0.6,42.3619
+0.61,41.8381
+0.62,40.2381
+0.63,39.1286
+0.64,38.7857
+0.65,37.081
+0.66,36.9524
+0.67,36.581
+0.68,36.281
+0.69,35.3476
+0.7,34.8905
+0.71,34.1667
+0.72,33.6714
+0.73,32.9667
+0.74,31.8429
+0.75,31.5429
+0.76,31.1476
+0.77,30.9905
+0.78,29.9571
+0.79,29.1333
+0.8,28.7857
+0.81,28.4429
+0.82,28.3476
+0.83,27.5429
+0.84,27.4333
+0.85,27.6048
+0.86,27.1762
+0.87,27.2
+0.88,26.4333
+0.89,25.7619
+0.9,24.8095
+0.91,24.7429
+0.92,24.2857
+0.93,24.1714
+0.94,23.5667
+0.95,23.5476
+0.96,23.3952
+0.97,22.919
+0.98,22.3095
+0.99,21.8048
+1,21.2857
+1.01,21.2048
+1.02,20.8429
+1.03,20.4429
+1.04,20.0048
+1.05,19.9381
+1.06,19.5
+1.07,19.8667
+1.08,18.9333
+1.09,19.1381
+1.1,18.9619
+1.11,18.5476
+1.12,17.9048
+1.13,17.7571
+1.14,18.5333
+1.15,18.3762
+1.16,18.3571
+1.17,18.3286
+1.18,18.2762
+1.19,18.3952
+1.2,17.5952
+1.21,18.1524
+1.22,18.1952
+1.23,17.8476
+1.24,17.9095
+1.25,17.5048
+1.26,17.5
+1.27,15.9619
+1.28,16.2095
+1.29,16.181
+1.3,15.6952
+1.31,15.7095
+1.32,15.4619
+1.33,15.9476
+1.34,16
+1.35,16.1952
+1.36,16.1143
+1.37,15.7429
+1.38,15.5762
+1.39,15.7048
+1.4,15.8095
+1.41,15.6667
+1.42,14.9048
+1.43,14.5857
+1.44,14.7524
+1.45,14.7571
+1.46,14.9762
+1.47,14.5333
+1.48,14.5524
+1.49,14.0143
+1.5,13.6286
+1.51,13.4429
+1.52,13.4667
+1.53,13.319
+1.54,12.9333
+1.55,13.1238
+1.56,12.7476
+1.57,12.9333
+1.58,13.0714
+1.59,13.0714
+1.6,12.7619
+1.61,12.4238
+1.62,12.5143
+1.63,12.9143
+1.64,12.5714
+1.65,13.3667
+1.66,13.2286
+1.67,13.7905
+1.68,13.7571
+1.69,13.5905
+1.7,12.9667
+1.71,12.981
+1.72,12.8857
+1.73,12.919
+1.74,13.0143
+1.75,13.0095
+1.76,12.3857
+1.77,12.5571
+1.78,12.3429
+1.79,12.7571
+1.8,12.681
+1.81,12.5429
+1.82,12.1857
+1.83,12.7905
+1.84,12.5571
+1.85,12.8429
+1.86,12.5476
+1.87,12.5714
+1.88,12.3762
+1.89,11.9952
+1.9,11.4571
+1.91,11.3
+1.92,11.1524
+1.93,11.681
+1.94,11.619
+1.95,11.9048
+1.96,12
+1.97,12.0762
+1.98,11.9143
+1.99,11.7619
+2,11.5333
diff --git a/test/kinetic_intensity_data.jl b/test/kinetic_intensity_data.jl
deleted file mode 100644
index 353015b..0000000
--- a/test/kinetic_intensity_data.jl
+++ /dev/null
@@ -1,202 +0,0 @@
-data = [
-    66.0952
-    104.762
-    110.333
-    114.905
-    122.238
-    125.429
-    125.429
-    123.476
-    121.286
-    118.857
-    117.667
-    116.143
-    113.857
-    111.571
-    108.81
-    105.952
-    104.048
-    102.048
-    100.143
-    98.5238
-    96.2381
-    94.381
-    91.6667
-    89.5714
-    87.1429
-    84.8571
-    83.4286
-    81.1905
-    78.9048
-    77.0476
-    75.4762
-    73.4762
-    71.8095
-    70.6667
-    68.381
-    67.3333
-    65.0952
-    63.7143
-    62.0476
-    60.8571
-    59.619
-    58.2857
-    57.4762
-    56.4762
-    55.8095
-    54.5238
-    53.0
-    51.8571
-    50.4286
-    49.381
-    47.9524
-    47.3714
-    46.8952
-    46.4857
-    45.9048
-    45.0762
-    44.3238
-    43.4143
-    43.5429
-    42.3619
-    41.8381
-    40.2381
-    39.1286
-    38.7857
-    37.081
-    36.9524
-    36.581
-    36.281
-    35.3476
-    34.8905
-    34.1667
-    33.6714
-    32.9667
-    31.8429
-    31.5429
-    31.1476
-    30.9905
-    29.9571
-    29.1333
-    28.7857
-    28.4429
-    28.3476
-    27.5429
-    27.4333
-    27.6048
-    27.1762
-    27.2
-    26.4333
-    25.7619
-    24.8095
-    24.7429
-    24.2857
-    24.1714
-    23.5667
-    23.5476
-    23.3952
-    22.919
-    22.3095
-    21.8048
-    21.2857
-    21.2048
-    20.8429
-    20.4429
-    20.0048
-    19.9381
-    19.5
-    19.8667
-    18.9333
-    19.1381
-    18.9619
-    18.5476
-    17.9048
-    17.7571
-    18.5333
-    18.3762
-    18.3571
-    18.3286
-    18.2762
-    18.3952
-    17.5952
-    18.1524
-    18.1952
-    17.8476
-    17.9095
-    17.5048
-    17.5
-    15.9619
-    16.2095
-    16.181
-    15.6952
-    15.7095
-    15.4619
-    15.9476
-    16.0
-    16.1952
-    16.1143
-    15.7429
-    15.5762
-    15.7048
-    15.8095
-    15.6667
-    14.9048
-    14.5857
-    14.7524
-    14.7571
-    14.9762
-    14.5333
-    14.5524
-    14.0143
-    13.6286
-    13.4429
-    13.4667
-    13.319
-    12.9333
-    13.1238
-    12.7476
-    12.9333
-    13.0714
-    13.0714
-    12.7619
-    12.4238
-    12.5143
-    12.9143
-    12.5714
-    13.3667
-    13.2286
-    13.7905
-    13.7571
-    13.5905
-    12.9667
-    12.981
-    12.8857
-    12.919
-    13.0143
-    13.0095
-    12.3857
-    12.5571
-    12.3429
-    12.7571
-    12.681
-    12.5429
-    12.1857
-    12.7905
-    12.5571
-    12.8429
-    12.5476
-    12.5714
-    12.3762
-    11.9952
-    11.4571
-    11.3
-    11.1524
-    11.681
-    11.619
-    11.9048
-    12.0
-    12.0762
-    11.9143
-    11.7619
-    11.5333
-]
\ No newline at end of file
diff --git a/test/runtests.jl b/test/runtests.jl
index 9c2b172..8e19776 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -1,3 +1,5 @@
+using CSV
+using DataFrames
 using EOptInterface
 using Ipopt
 using JuMP
@@ -191,12 +193,11 @@ end
 
     @mtkcompile system = KineticParameterEstimation()
 
-    tspan = (0.0, 2.0)
-    tstep = 0.01
-    N = Int(floor((tspan[2] - tspan[1])/tstep)) + 1
-
-    include("kinetic_intensity_data.jl")
+    data = CSV.read("kinetic_intensity_data.csv", DataFrame)
     intensity(x_A, x_B, x_D) = x_A + 2/21*x_B + 2/21*x_D
+    tspan = (data.time[1], data.time[end])
+    tstep = data.time[2] - data.time[1]
+    N = length(data.time)
 
     model = JuMP.Model(Ipopt.Optimizer)
     V = length(unknowns(system))
@@ -209,7 +210,7 @@ end
 
     EOptInterface.register_odesystem(model, system, tspan, tstep, "EE")
 
-    JuMP.@objective(model, Min, sum((intensity(z[5,i], z[4,i], z[3,i]) - data[i-1])^2 for i in 2:N))
+    JuMP.@objective(model, Min, sum((intensity(z[5,i], z[4,i], z[3,i]) - data.intensity[i])^2 for i in 1:N))
     JuMP.optimize!(model)
 
     @test JuMP.termination_status(model) == JuMP.LOCALLY_SOLVED
@@ -223,7 +224,7 @@ end
     pU = [1200.0, 1200.0, 40.0]
     JuMP.@variable(model, pL[i] <= p[i=1:3] <= pU[i])
     EOptInterface.register_odesystem(model, system, tspan, tstep, "IE")
-    JuMP.@objective(model, Min, sum((intensity(z[5,i], z[4,i], z[3,i]) - data[i-1])^2 for i in 2:N))
+    JuMP.@objective(model, Min, sum((intensity(z[5,i], z[4,i], z[3,i]) - data.intensity[i])^2 for i in 1:N))
     JuMP.optimize!(model)
 
     @test JuMP.termination_status(model) == JuMP.LOCALLY_SOLVED