diff --git a/Project.toml b/Project.toml
index 6ff0a08..b4c7d9e 100644
--- a/Project.toml
+++ b/Project.toml
@@ -4,6 +4,7 @@ authors = ["Joseph Choi <jsphchoi@mit.edu>", "Dimitri Alston <dimitri.alston@uco
 version = "0.2.0"
 
 [deps]
+DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae"
 JuMP = "4076af6c-e467-56ae-b986-b466b2749572"
 ModelingToolkit = "961ee093-0014-501f-94e3-6117800e7a78"
 Reexport = "189a3867-3050-52da-a836-e630ba90ab69"
@@ -12,6 +13,7 @@ SymbolicUtils = "d1185830-fcd6-423d-90d6-eec64667417b"
 Symbolics = "0c5d862f-8b57-4792-8d23-62f2024744c7"
 
 [compat]
+DocStringExtensions = "0.9"
 JuMP = "1"
 ModelingToolkit = "11"
 Reexport = "1"
diff --git a/README.md b/README.md
index c8a18cd..07a64f1 100644
--- a/README.md
+++ b/README.md
@@ -10,32 +10,9 @@
 |:-----------------------------------------------------------------:|
 | [![](https://img.shields.io/badge/docs-stable-blue.svg)](https://PSORLab.github.io/EOptInterface.jl/stable) [![](https://img.shields.io/badge/docs-latest-blue.svg)](https://PSORLab.github.io/EOptInterface.jl/dev) |
 
-## Feature Summary
-
-```julia
-decision_vars(::ModelingToolkit.System)
-```
-Returns the decision variables for an optimization problem formulated from a `ModelingToolkit` system.
-
-```julia
-register_nlsystem(::JuMP.Model, ::ModelingToolkit.System, obj::Symbolics.Num, ineqs::Vector{Symbolics.Num})
-```
-Automatically formulates algebraic `JuMP` constraints and objective function from
-an algebraic `ModelingToolkit` system and user-provided constraints and objective symbolic expressions.
-
-```julia
-full_solution(::JuMP.Model, ::ModelingToolkit.System)
-```
-Returns a dictionary of optimal solution values for the observed variables of an algebraic `ModelingToolkit` system if the `JuMP` model is solved.
-
-```julia
-register_odesystem(::JuMP.Model, ::ModelingToolkit.System, tspan::Tuple{Real,Real}, tstep::Real, solver::String)
-```
-Automatically applies forward transcription and registers the discretized ODE `ModelingToolkit` system as algebraic `JuMP` constraints. Available integration methods: `"EE"` (explicit Euler), `"IE"` (implicit Euler).
-
 ## Examples
 
-The code for these examples can be found in the [`examples/`](https://github.com/PSORLab/EOptInterface.jl/tree/main/examples) subdirectory.
+The code and notebooks for these examples can be found in the [`examples/`](https://github.com/PSORLab/EOptInterface.jl/tree/main/examples) subdirectory.
 
 ### [Algebraic System](https://github.com/PSORLab/EOptInterface.jl/blob/main/examples/algebraic_model.jl)
 
@@ -43,7 +20,7 @@ An optimal reactor-separator-recycle process design problem originally presented
 
 ### [ODE System](https://github.com/PSORLab/EOptInterface.jl/blob/main/examples/ode_model.jl)
 
-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` [[6](#references)].
+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)
diff --git a/src/EOptInterface.jl b/src/EOptInterface.jl
index 3669e88..b9dcf1c 100644
--- a/src/EOptInterface.jl
+++ b/src/EOptInterface.jl
@@ -13,12 +13,13 @@
 module EOptInterface
 
     # Imports
+    import DocStringExtensions
     import JuMP
     import ModelingToolkit
     import ModelingToolkit: t_nounits as t, D_nounits as D
     import Reexport
-    import Symbolics
     import SymbolicUtils
+    import Symbolics
 
     # Reexports
     Reexport.@reexport using SciCompDSL
diff --git a/src/userfuncs.jl b/src/userfuncs.jl
index 26497ae..fa237f6 100644
--- a/src/userfuncs.jl
+++ b/src/userfuncs.jl
@@ -13,9 +13,10 @@
 ################################################################################
 
 """
-    decision_vars(sys)
+    $(DocStringExtensions.TYPEDSIGNATURES)
 
-Returns the decision variables for an optimization problem formulated from a ModelingToolkit system.
+Returns the decision variables for an optimization problem from a 
+`ModelingToolkit.System`.
 """
 function decision_vars(sys::ModelingToolkit.System)
     return [
@@ -25,16 +26,11 @@ function decision_vars(sys::ModelingToolkit.System)
 end
 
 """
-    register_nlsystem(model, sys, obj, ineqs)
+    $(DocStringExtensions.TYPEDSIGNATURES)
 
-Automatically formulates algebraic JuMP constraints and objective function from
-an algebraic ModelingToolkit system and user-provided constraints and objective symbolic expressions.
-
-# Arguments
-- `model::JuMP.Model`: the JuMP model
-- `sys::ModelingToolkit.System`: the ModelingToolkit system
-- `obj::Symbolics.Num`: a symbolic expression of the objective function using the ModelingToolkit system variables
-- `ineqs::Vector{Symbolics.Num}`: a vector of symbolic expressions of inequality constraints using the ModelingToolkit system variables
+Automatically formulates and adds user-provided `Symbolics.Num` objective 
+function and `Vector{Symbolics.Num}` constraints from a `ModelingToolkit.System` 
+to a `JuMP.Model`.
 """
 function register_nlsystem(model::JuMP.Model, sys::ModelingToolkit.System, obj::Symbolics.Num, ineqs::Vector{Symbolics.Num})
     h = EOptInterface.mtk_generate_model_equations(sys)
@@ -51,49 +47,44 @@ function register_nlsystem(model::JuMP.Model, sys::ModelingToolkit.System, obj::
 end
 
 """
-    register_odesystem(model, sys, tspan, tstep, integrator)
-
-Automatically applies forward transcription and registers the discretized ODE ModelingToolkit system as algebraic JuMP constraints.
+    $(DocStringExtensions.TYPEDSIGNATURES)
 
-# Arguments
-- `model::JuMP.Model`: the JuMP model
-- `sys::ModelingToolkit.System`: the ModelingToolkit system
-- `tspan::Tuple{Real,Real}`: the time span over which the dynamic model is simulated
-- `tstep::Real`: the time step used in the integration scheme
-- `integrator::String`: integration scheme used in discretization, `"EE"` for explicit Euler or `"IE"` for implicit Euler
+Automatically applies specified direct transcription method and registers the 
+discretized ODE `ModelingToolkit.System` as algebraic JuMP constraints.
+Current supports Explicit Euler "EE" and Implicit Euler "IE".
 """
-function register_odesystem(model::JuMP.Model, odesys::ModelingToolkit.System, tspan::Tuple{Real,Real}, tstep::Real, integrator::String)
+function register_odesystem(model::JuMP.Model, sys::ModelingToolkit.System, tspan::Tuple{Real,Real}, tstep::Real, integrator::String)
     if integrator != "EE" && integrator != "IE"
         error("Available integrators: EE, IE")
     end
     # Number of discrete time nodes
     N = Int(floor((tspan[2] - tspan[1])/tstep)) + 1
     # Number of ODE variables
-    V = length(ModelingToolkit.unknowns(odesys))
-    param_dict = copy(ModelingToolkit.initial_conditions(odesys).dict)
-    for var in ModelingToolkit.unknowns(odesys)
+    V = length(ModelingToolkit.unknowns(sys))
+    param_dict = copy(ModelingToolkit.initial_conditions(sys).dict)
+    for var in ModelingToolkit.unknowns(sys)
         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))
+        dxj_expr = ModelingToolkit.full_equations(sys)[j].rhs
+        dxj_expr = SymbolicUtils.substitute(dxj_expr, ModelingToolkit.bindings(sys))
+        dxj_expr = SymbolicUtils.substitute(dxj_expr, ModelingToolkit.bindings(sys))
         # 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)
         end
         dxj = Symbolics.build_function(
             dxj_expr,
-            EOptInterface.decision_vars(odesys)..., 
+            EOptInterface.decision_vars(sys)..., 
             expression = Val{false}
             )
         push!(dx, dxj)
     end
-    ps = JuMP.all_variables(model)[end-length(setdiff(EOptInterface.decision_vars(odesys), ModelingToolkit.unknowns(odesys)))+1:end]
+    ps = JuMP.all_variables(model)[end-length(setdiff(EOptInterface.decision_vars(sys), ModelingToolkit.unknowns(sys)))+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.initial_conditions(odesys)[ModelingToolkit.unknowns(odesys)[i]].val for i in eachindex(ModelingToolkit.unknowns(odesys))], force=true)
+    JuMP.fix.(xs[:,1], [ModelingToolkit.initial_conditions(sys)[ModelingToolkit.unknowns(sys)[i]].val for i in eachindex(ModelingToolkit.unknowns(sys))], 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...))
@@ -104,9 +95,10 @@ function register_odesystem(model::JuMP.Model, odesys::ModelingToolkit.System, t
 end
 
 """
-    full_solution(model, sys)
+    $(DocStringExtensions.TYPEDSIGNATURES)
 
-Returns a dictionary of optimal solution values for the observed variables of an algebraic ModelingToolkit system if the JuMP model is solved.
+Returns a dictionary of optimal solution values for the observed variables of a 
+`ModelingToolkit.System` if the `JuMP.Model` is solved.
 """
 function full_solution(model::JuMP.Model, sys::ModelingToolkit.System)
     vars = EOptInterface.decision_vars(sys)