diff --git a/src/ModelingToolkit.jl b/src/ModelingToolkit.jl
index 5b13325c62..9276ade371 100644
--- a/src/ModelingToolkit.jl
+++ b/src/ModelingToolkit.jl
@@ -188,7 +188,6 @@ include("modelingtoolkitize/nonlinearproblem.jl")
 
 include("systems/nonlinear/homotopy_continuation.jl")
 include("systems/nonlinear/initializesystem.jl")
-include("systems/diffeqs/first_order_transform.jl")
 include("systems/diffeqs/basic_transformations.jl")
 
 include("systems/pde/pdesystem.jl")
@@ -291,9 +290,8 @@ export isinput, isoutput, getbounds, hasbounds, getguess, hasguess, isdisturbanc
        tunable_parameters, isirreducible, getdescription, hasdescription,
        hasunit, getunit, hasconnect, getconnect,
        hasmisc, getmisc, state_priority
-export ode_order_lowering, dae_order_lowering, liouville_transform,
-       change_independent_variable, substitute_component, add_accumulations,
-       noise_to_brownians
+export liouville_transform, change_independent_variable, substitute_component,
+       add_accumulations, noise_to_brownians
 export PDESystem
 export Differential, expand_derivatives, @derivatives
 export Equation, ConstrainedEquation
diff --git a/src/systems/diffeqs/first_order_transform.jl b/src/systems/diffeqs/first_order_transform.jl
deleted file mode 100644
index d017ea362b..0000000000
--- a/src/systems/diffeqs/first_order_transform.jl
+++ /dev/null
@@ -1,106 +0,0 @@
-"""
-$(TYPEDSIGNATURES)
-
-Takes a Nth order System and returns a new System written in first order
-form by defining new variables which represent the N-1 derivatives.
-"""
-function ode_order_lowering(sys::System)
-    iv = get_iv(sys)
-    eqs_lowered, new_vars = ode_order_lowering(equations(sys), iv, unknowns(sys))
-    @set! sys.eqs = eqs_lowered
-    @set! sys.unknowns = new_vars
-    return sys
-end
-
-function dae_order_lowering(sys::System)
-    iv = get_iv(sys)
-    eqs_lowered, new_vars = dae_order_lowering(equations(sys), iv, unknowns(sys))
-    @set! sys.eqs = eqs_lowered
-    @set! sys.unknowns = new_vars
-    return sys
-end
-
-function ode_order_lowering(eqs, iv, unknown_vars)
-    var_order = OrderedDict{Any, Int}()
-    D = Differential(iv)
-    diff_eqs = Equation[]
-    diff_vars = []
-    alge_eqs = Equation[]
-
-    for (i, eq) in enumerate(eqs)
-        if !isdiffeq(eq)
-            push!(alge_eqs, eq)
-        else
-            var, maxorder = var_from_nested_derivative(eq.lhs)
-            maxorder > get(var_order, var, 1) && (var_order[var] = maxorder)
-            var′ = lower_varname(var, iv, maxorder - 1)
-            rhs′ = diff2term_with_unit(eq.rhs, iv)
-            push!(diff_vars, var′)
-            push!(diff_eqs, D(var′) ~ rhs′)
-        end
-    end
-
-    for (var, order) in var_order
-        for o in (order - 1):-1:1
-            lvar = lower_varname(var, iv, o - 1)
-            rvar = lower_varname(var, iv, o)
-            push!(diff_vars, lvar)
-
-            rhs = rvar
-            eq = Differential(iv)(lvar) ~ rhs
-            push!(diff_eqs, eq)
-        end
-    end
-
-    # we want to order the equations and variables to be `(diff, alge)`
-    return (vcat(diff_eqs, alge_eqs), vcat(diff_vars, setdiff(unknown_vars, diff_vars)))
-end
-
-function dae_order_lowering(eqs, iv, unknown_vars)
-    var_order = OrderedDict{Any, Int}()
-    D = Differential(iv)
-    diff_eqs = Equation[]
-    diff_vars = OrderedSet()
-    alge_eqs = Equation[]
-    vars = Set()
-    subs = Dict()
-
-    for (i, eq) in enumerate(eqs)
-        vars!(vars, eq)
-        n_diffvars = 0
-        for vv in vars
-            isdifferential(vv) || continue
-            var, maxorder = var_from_nested_derivative(vv)
-            isparameter(var) && continue
-            n_diffvars += 1
-            order = get(var_order, var, nothing)
-            seen = order !== nothing
-            if !seen
-                order = 1
-            end
-            maxorder > order && (var_order[var] = maxorder)
-            var′ = lower_varname(var, iv, maxorder - 1)
-            subs[vv] = D(var′)
-            if !seen
-                push!(diff_vars, var′)
-            end
-        end
-        n_diffvars == 0 && push!(alge_eqs, eq)
-        empty!(vars)
-    end
-
-    for (var, order) in var_order
-        for o in (order - 1):-1:1
-            lvar = lower_varname(var, iv, o - 1)
-            rvar = lower_varname(var, iv, o)
-            push!(diff_vars, lvar)
-
-            rhs = rvar
-            eq = Differential(iv)(lvar) ~ rhs
-            push!(diff_eqs, eq)
-        end
-    end
-
-    return ([diff_eqs; substitute.(eqs, (subs,))],
-        vcat(collect(diff_vars), setdiff(unknown_vars, diff_vars)))
-end
diff --git a/test/lowering_solving.jl b/test/lowering_solving.jl
deleted file mode 100644
index 280848f6ad..0000000000
--- a/test/lowering_solving.jl
+++ /dev/null
@@ -1,76 +0,0 @@
-using ModelingToolkit, OrdinaryDiffEq, Test, LinearAlgebra
-using ModelingToolkit: t_nounits as t, D_nounits as D
-
-@parameters σ ρ β
-@variables x(t) y(t) z(t) k(t)
-
-eqs = [D(D(x)) ~ σ * (y - x),
-    D(y) ~ x * (ρ - z) - y,
-    D(z) ~ x * y - β * z]
-
-@named sys′ = System(eqs, t)
-sys = ode_order_lowering(sys′)
-
-eqs2 = [0 ~ x * y - k,
-    D(D(x)) ~ σ * (y - x),
-    D(y) ~ x * (ρ - z) - y,
-    D(z) ~ x * y - β * z]
-@named sys2 = System(eqs2, t, [x, y, z, k], parameters(sys′))
-sys2 = ode_order_lowering(sys2)
-# test equation/variable ordering
-ModelingToolkit.calculate_massmatrix(sys2) == Diagonal([1, 1, 1, 1, 0])
-
-u0 = [D(x) => 2.0,
-    x => 1.0,
-    y => 0.0,
-    z => 0.0]
-
-p = [σ => 28.0,
-    ρ => 10.0,
-    β => 8 / 3]
-
-tspan = (0.0, 100.0)
-
-sys = complete(sys)
-prob = ODEProblem(sys, u0, tspan, p, jac = true)
-probexpr = ODEProblem(sys, u0, tspan, p; jac = true, expression = Val{true})
-sol = solve(prob, Tsit5())
-solexpr = solve(eval(prob), Tsit5())
-@test all(x -> x == 0, Array(sol - solexpr))
-#using Plots; plot(sol,idxs=(:x,:y))
-
-@parameters σ ρ β
-@variables x(t) y(t) z(t)
-
-eqs = [D(x) ~ σ * (y - x),
-    D(y) ~ x * (ρ - z) - y,
-    D(z) ~ x * y - β * z]
-
-lorenz1 = System(eqs, t, name = :lorenz1)
-lorenz2 = System(eqs, t, name = :lorenz2)
-
-@variables α(t)
-@parameters γ
-connections = [0 ~ lorenz1.x + lorenz2.y + α * γ]
-@named connected = System(connections, t, [α], [γ], systems = [lorenz1, lorenz2])
-connected = complete(connected)
-u0 = [lorenz1.x => 1.0,
-    lorenz1.y => 0.0,
-    lorenz1.z => 0.0,
-    lorenz2.x => 0.0,
-    lorenz2.y => 1.0,
-    lorenz2.z => 0.0]
-
-p = [lorenz1.σ => 10.0,
-    lorenz1.ρ => 28.0,
-    lorenz1.β => 8 / 3,
-    lorenz2.σ => 10.0,
-    lorenz2.ρ => 28.0,
-    lorenz2.β => 8 / 3,
-    γ => 2.0]
-
-tspan = (0.0, 100.0)
-prob = ODEProblem(connected, u0, tspan, p)
-sol = solve(prob, Rodas5())
-@test maximum(sol[lorenz1.x] + sol[lorenz2.y] + 2sol[α]) < 1e-12
-#using Plots; plot(sol,idxs=(:α,Symbol(lorenz1.x),Symbol(lorenz2.y)))
diff --git a/test/odesystem.jl b/test/odesystem.jl
index 9c8a67e7eb..2bfca7dcdd 100644
--- a/test/odesystem.jl
+++ b/test/odesystem.jl
@@ -153,32 +153,6 @@ du = [0.0]
 f(du, [1.0], [t -> t + 2], 5.0)
 @test du ≈ [27561]
 
-@testset "Issue#17: Conversion to first order ODEs" begin
-    D3 = D^3
-    D2 = D^2
-    @variables u(t) uˍtt(t) uˍt(t) xˍt(t)
-    eqs = [D3(u) ~ 2(D2(u)) + D(u) + D(x) + 1
-           D2(x) ~ D(x) + 2]
-    @named de = System(eqs, t)
-    de1 = ode_order_lowering(de)
-
-    @testset "Issue#219: Ordering of equations in `ode_order_lowering`" begin
-        lowered_eqs = [D(uˍtt) ~ 2uˍtt + uˍt + xˍt + 1
-                       D(xˍt) ~ xˍt + 2
-                       D(uˍt) ~ uˍtt
-                       D(u) ~ uˍt
-                       D(x) ~ xˍt]
-        @test isequal(
-            [ModelingToolkit.var_from_nested_derivative(eq.lhs)[1] for eq in equations(de1)],
-            unknowns(@named lowered = System(lowered_eqs, t)))
-    end
-
-    test_diffeq_inference("first-order transform", de1, t, [uˍtt, xˍt, uˍt, u, x], [])
-    du = zeros(5)
-    ODEFunction(complete(de1))(du, ones(5), nothing, 0.1)
-    @test du == [5.0, 3.0, 1.0, 1.0, 1.0]
-end
-
 # Internal calculations
 @parameters σ
 a = y - x
@@ -348,16 +322,6 @@ eqs = [D(x) ~ σ * (y - x),
     @test issym(equations(sys)[1].rhs)
 end
 
-@testset "Issue#708" begin
-    @parameters a
-    @variables x(t) y(t) z(t)
-    @named sys = System([D(x) ~ y, 0 ~ x + z, 0 ~ x - y], t, [z, y, x], [])
-
-    sys2 = ode_order_lowering(sys)
-    M = ModelingToolkit.calculate_massmatrix(sys2)
-    @test M == Diagonal([1, 0, 0])
-end
-
 # issue #609
 @variables x1(t) x2(t)
 
@@ -416,22 +380,6 @@ eqs = [
 ]
 @test_throws ArgumentError ModelingToolkit.System(eqs, t, vars, pars, name = :foo)
 
-@variables x(t)
-@parameters M b k
-eqs = [D(D(x)) ~ -b / M * D(x) - k / M * x]
-ps = [M, b, k]
-default_u0 = [D(x) => 0.0, x => 10.0]
-default_p = [M => 1.0, b => 1.0, k => 1.0]
-@named sys = System(eqs, t, [x], ps; defaults = [default_u0; default_p])
-sys = ode_order_lowering(sys)
-sys = complete(sys)
-prob = ODEProblem(sys, nothing, tspan)
-sol = solve(prob, Tsit5())
-@test sol.t[end] == tspan[end]
-@test sum(abs, sol.u[end]) < 1
-prob = ODEProblem{false}(sys, nothing, tspan; u0_constructor = x -> SVector(x...))
-@test prob.u0 isa SVector
-
 # check_eqs_u0 kwarg test
 @variables x1(t) x2(t)
 eqs = [D(x1) ~ -x1]
@@ -1530,64 +1478,6 @@ end
     @test osys1 !== osys2
 end
 
-@testset "dae_order_lowering basic test" begin
-    @parameters a
-    @variables x(t) y(t) z(t)
-    @named dae_sys = System([
-            D(x) ~ y,
-            0 ~ x + z,
-            0 ~ x - y + z
-        ], t, [z, y, x], [])
-
-    lowered_dae_sys = dae_order_lowering(dae_sys)
-    @variables x1(t) y1(t) z1(t)
-    expected_eqs = [
-        0 ~ x + z,
-        0 ~ x - y + z,
-        Differential(t)(x) ~ y
-    ]
-    lowered_eqs = equations(lowered_dae_sys)
-    sorted_lowered_eqs = sort(lowered_eqs, by = string)
-    sorted_expected_eqs = sort(expected_eqs, by = string)
-    @test sorted_lowered_eqs == sorted_expected_eqs
-
-    expected_vars = Set([z, y, x])
-    lowered_vars = Set(unknowns(lowered_dae_sys))
-    @test lowered_vars == expected_vars
-end
-
-@testset "dae_order_lowering test with structural_simplify" begin
-    @variables x(t) y(t) z(t)
-    @parameters M b k
-    eqs = [
-        D(D(x)) ~ -b / M * D(x) - k / M * x,
-        0 ~ y - D(x),
-        0 ~ z - x
-    ]
-    ps = [M, b, k]
-    default_u0 = [
-        D(x) => 0.0, x => 10.0, y => 0.0, z => 10.0
-    ]
-    default_p = [M => 1.0, b => 1.0, k => 1.0]
-    @named dae_sys = System(eqs, t, [x, y, z], ps; defaults = [default_u0; default_p])
-
-    simplified_dae_sys = structural_simplify(dae_sys)
-
-    lowered_dae_sys = dae_order_lowering(simplified_dae_sys)
-    lowered_dae_sys = complete(lowered_dae_sys)
-
-    tspan = (0.0, 10.0)
-    prob = ODEProblem(lowered_dae_sys, nothing, tspan)
-    sol = solve(prob, Tsit5())
-
-    @test sol.t[end] == tspan[end]
-    @test sum(abs, sol.u[end]) < 1
-
-    prob = ODEProblem{false}(
-        lowered_dae_sys, nothing, tspan; u0_constructor = x -> SVector(x...))
-    @test prob.u0 isa SVector
-end
-
 @testset "Constraint system construction" begin
     @variables x(..) y(..) z(..)
     @parameters a b c d e
diff --git a/test/runtests.jl b/test/runtests.jl
index 1b8b5e03db..9a184df054 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -51,7 +51,6 @@ end
             @safetestset "Symbolic Event Test" include("symbolic_events.jl")
             @safetestset "Stream Connect Test" include("stream_connectors.jl")
             @safetestset "Domain Connect Test" include("domain_connectors.jl")
-            @safetestset "Lowering Integration Test" include("lowering_solving.jl")
             @safetestset "Dependency Graph Test" include("dep_graphs.jl")
             @safetestset "Function Registration Test" include("function_registration.jl")
             @safetestset "Precompiled Modules Test" include("precompile_test.jl")
diff --git a/test/state_selection.jl b/test/state_selection.jl
index 6db8e8c5a0..ba74cc04a2 100644
--- a/test/state_selection.jl
+++ b/test/state_selection.jl
@@ -23,7 +23,6 @@ end
 @test_skip let pss = partial_state_selection(sys)
     @test length(equations(pss)) == 1
     @test length(unknowns(pss)) == 2
-    @test length(equations(ode_order_lowering(pss))) == 2
 end
 
 @parameters σ ρ β
diff --git a/test/structural_transformation/index_reduction.jl b/test/structural_transformation/index_reduction.jl
index 26f73db4e7..39c989aa3f 100644
--- a/test/structural_transformation/index_reduction.jl
+++ b/test/structural_transformation/index_reduction.jl
@@ -3,27 +3,10 @@ using Graphs
 using DiffEqBase
 using Test
 using UnPack
+using OrdinaryDiffEq
+using LinearAlgebra
 using ModelingToolkit: t_nounits as t, D_nounits as D
 
-# Define some variables
-@parameters L g
-@variables x(t) y(t) w(t) z(t) T(t) xˍt(t) yˍt(t) xˍˍt(t) yˍˍt(t)
-
-eqs2 = [D(D(x)) ~ T * x,
-    D(D(y)) ~ T * y - g,
-    0 ~ x^2 + y^2 - L^2]
-pendulum2 = System(eqs2, t, [x, y, T], [L, g], name = :pendulum)
-lowered_sys = ModelingToolkit.ode_order_lowering(pendulum2)
-
-lowered_eqs = [D(xˍt) ~ T * x,
-    D(yˍt) ~ T * y - g,
-    D(x) ~ xˍt,
-    D(y) ~ yˍt,
-    0 ~ x^2 + y^2 - L^2]
-@test System(lowered_eqs, t, [xˍt, yˍt, x, y, T], [L, g], name = :pendulum) ==
-      lowered_sys
-@test isequal(equations(lowered_sys), lowered_eqs)
-
 # Simple pendulum in cartesian coordinates
 eqs = [D(x) ~ w,
     D(y) ~ z,
@@ -39,50 +22,6 @@ state = TearingState(pendulum)
       map(x -> x == 0 ? StructuralTransformations.unassigned : x,
     [3, 4, 2, 5, 0, 0, 0, 0, 0])
 
-using ModelingToolkit
-@parameters L g
-@variables x(t) y(t) w(t) z(t) T(t) xˍt(t) yˍt(t)
-idx1_pendulum = [D(x) ~ w,
-    D(y) ~ z,
-    #0 ~ x^2 + y^2 - L^2,
-    D(w) ~ T * x,
-    D(z) ~ T * y - g,
-    # intermediate 1: 0 ~ 2x*D(x) + 2y*D(y) - 0,
-    # intermediate 2(a): 0 ~ 2x*w + 2y*z - 0, (substitute D(x) and D(y))
-    #0 ~ 2x*w + 2y*z,
-    # D(D(x)) ~ D(w) and substitute the rhs
-    D(xˍt) ~ T * x,
-    # D(D(y)) ~ D(z) and substitute the rhs
-    D(yˍt) ~ T * y - g,
-    # 2x*D(D(x)) + 2*D(x)*D(x) + 2y*D(D(y)) + 2*D(y)*D(y) and
-    # substitute the rhs
-    0 ~ 2x * (T * x) + 2 * xˍt * xˍt + 2y * (T * y - g) + 2 * yˍt * yˍt]
-@named idx1_pendulum = System(idx1_pendulum, t, [x, y, w, z, xˍt, yˍt, T], [L, g])
-first_order_idx1_pendulum = complete(ode_order_lowering(idx1_pendulum))
-
-using OrdinaryDiffEq
-using LinearAlgebra
-prob = ODEProblem(first_order_idx1_pendulum,
-    #  [x, y, w, z, xˍt, yˍt, T]
-    [1, 0, 0, 0, 0, 0, 0.0],# 0, 0, 0, 0],
-    (0, 10.0),
-    [1, 9.8])
-sol = solve(prob, Rodas5());
-#plot(sol, idxs=(1, 2))
-
-new_sys = complete(dae_index_lowering(ModelingToolkit.ode_order_lowering(pendulum2)))
-
-prob_auto = ODEProblem(new_sys,
-    [D(x) => 0,
-        D(y) => 0,
-        x => 1,
-        y => 0,
-        T => 0.0],
-    (0, 100.0),
-    [L => 1, g => 9.8])
-sol = solve(prob_auto, Rodas5());
-#plot(sol, idxs=(x, y))
-
 # Define some variables
 @parameters L g
 @variables x(t) y(t) T(t)
@@ -92,29 +31,6 @@ eqs2 = [D(D(x)) ~ T * x,
     0 ~ x^2 + y^2 - L^2]
 pendulum2 = System(eqs2, t, [x, y, T], [L, g], name = :pendulum)
 
-# Turn into a first order differential equation system
-first_order_sys = ModelingToolkit.ode_order_lowering(pendulum2)
-
-# Perform index reduction to get an Index 1 DAE
-new_sys = complete(dae_index_lowering(first_order_sys))
-
-u0 = [
-    D(x) => 0.0,
-    D(y) => 0.0,
-    x => 1.0,
-    y => 0.0,
-    T => 0.0
-]
-
-p = [
-    L => 1.0,
-    g => 9.8
-]
-
-prob_auto = ODEProblem(new_sys, u0, (0.0, 10.0), p)
-sol = solve(prob_auto, Rodas5());
-#plot(sol, idxs=(D(x), y))
-
 @test_skip begin
     let pss_pendulum2 = partial_state_selection(pendulum2)
         length(equations(pss_pendulum2)) <= 6