diff --git a/src/ModelingToolkit.jl b/src/ModelingToolkit.jl
index 7a6b6e16a2..46fa3f3299 100644
--- a/src/ModelingToolkit.jl
+++ b/src/ModelingToolkit.jl
@@ -4,6 +4,7 @@ export Operation, Expression
 export calculate_jacobian, generate_jacobian, generate_function
 export independent_variables, dependent_variables, parameters
 export @register
+export modelingtoolkitize
 
 
 using DiffEqBase
diff --git a/src/utils.jl b/src/utils.jl
index 75f729c3f1..6aff728232 100644
--- a/src/utils.jl
+++ b/src/utils.jl
@@ -84,3 +84,27 @@ function vars!(vars, O)
 
     return vars
 end
+
+"""
+$(SIGNATURES)
+
+Generate `ODESystem`, dependent variables, and parameters from an `ODEProblem`.
+"""
+function modelingtoolkitize(prob::DiffEqBase.ODEProblem)
+    t, = @parameters t; vars = [Variable(Symbol(:x, i))(t) for i in eachindex(prob.u0)]; params = [Variable(Symbol(:α, i); known = true)() for i in eachindex(prob.p)];
+    D, = @derivatives D'~t
+
+    rhs = [D(var) for var in vars]
+
+    if DiffEqBase.isinplace(prob)
+        lhs = similar(vars, Any)
+        prob.f(lhs, vars, params, t)
+    else
+        lhs = prob.f(vars, params, t)
+    end
+
+    eqs = vcat([rhs[i] ~ lhs[i] for i in eachindex(prob.u0)]...)
+    de = ODESystem(eqs)
+
+    de, vars, params
+end
diff --git a/test/system_construction.jl b/test/system_construction.jl
index 154735b415..f11e579896 100644
--- a/test/system_construction.jl
+++ b/test/system_construction.jl
@@ -1,4 +1,5 @@
 using ModelingToolkit, StaticArrays, LinearAlgebra
+using DiffEqBase
 using Test
 
 # Define some variables
@@ -197,3 +198,26 @@ ns = NonlinearSystem(eqs, [x,y,z])
 nlsys_func = generate_function(ns, [x,y,z], [σ,ρ,β])
 jac = calculate_jacobian(ns)
 jac = generate_jacobian(ns)
+
+function lotka(u,p,t)
+  x = u[1]
+  y = u[2]
+  [p[1]*x - p[2]*x*y,
+  -p[3]*y + p[4]*x*y]
+end
+
+prob = ODEProblem(ODEFunction{false}(lotka),[1.0,1.0],(0.0,1.0),[1.5,1.0,3.0,1.0])
+de, vars, params = modelingtoolkitize(prob)
+ODEFunction(de, vars, params)(similar(prob.u0), prob.u0, prob.p, 0.1)
+
+function lotka(du,u,p,t)
+  x = u[1]
+  y = u[2]
+  du[1] = p[1]*x - p[2]*x*y
+  du[2] = -p[3]*y + p[4]*x*y
+end
+
+prob = ODEProblem(lotka,[1.0,1.0],(0.0,1.0),[1.5,1.0,3.0,1.0])
+
+de, vars, params = modelingtoolkitize(prob)
+ODEFunction(de, vars, params)(similar(prob.u0), prob.u0, prob.p, 0.1)