diff --git a/src/ModelingToolkit.jl b/src/ModelingToolkit.jl
index b0197245f9..50da559aaa 100644
--- a/src/ModelingToolkit.jl
+++ b/src/ModelingToolkit.jl
@@ -151,7 +151,7 @@ include("systems/optimization/modelingtoolkitize.jl")
 include("systems/pde/pdesystem.jl")
 
 include("systems/sparsematrixclil.jl")
-include("systems/discrete_system/discrete_system.jl")
+
 include("systems/unit_check.jl")
 include("systems/validation.jl")
 include("systems/dependency_graphs.jl")
@@ -198,7 +198,7 @@ export NonlinearProblem, BlockNonlinearProblem, NonlinearProblemExpr
 export OptimizationProblem, OptimizationProblemExpr, constraints
 export AutoModelingToolkit
 export SteadyStateProblem, SteadyStateProblemExpr
-export JumpProblem, DiscreteProblem
+export JumpProblem
 export NonlinearSystem, OptimizationSystem, ConstraintsSystem
 export alias_elimination, flatten
 export connect, domain_connect, @connector, Connection, Flow, Stream, instream
@@ -216,9 +216,6 @@ export SymScope, LocalScope, ParentScope, DelayParentScope, GlobalScope
 export independent_variable, equations, controls,
     observed, structure, full_equations
 export structural_simplify, expand_connections, linearize, linearization_function
-export DiscreteSystem,
-    DiscreteProblem, DiscreteProblemExpr, DiscreteFunction,
-    DiscreteFunctionExpr
 
 export calculate_jacobian, generate_jacobian, generate_function
 export calculate_control_jacobian, generate_control_jacobian
@@ -229,7 +226,6 @@ export calculate_hessian, generate_hessian
 export calculate_massmatrix, generate_diffusion_function
 export stochastic_integral_transform
 export TearingState, StateSelectionState
-export generate_difference_cb
 
 export BipartiteGraph, equation_dependencies, variable_dependencies
 export eqeq_dependencies, varvar_dependencies
diff --git a/src/clock.jl b/src/clock.jl
index 292577ee49..0ce64980f1 100644
--- a/src/clock.jl
+++ b/src/clock.jl
@@ -42,7 +42,7 @@ end
 has_time_domain(x::Num) = has_time_domain(value(x))
 has_time_domain(x) = false
 
-for op in [Differential, Difference]
+for op in [Differential]
     @eval input_timedomain(::$op, arg = nothing) = Continuous()
     @eval output_timedomain(::$op, arg = nothing) = Continuous()
 end
diff --git a/src/structural_transformation/codegen.jl b/src/structural_transformation/codegen.jl
index 1123420a87..1e154e584d 100644
--- a/src/structural_transformation/codegen.jl
+++ b/src/structural_transformation/codegen.jl
@@ -543,8 +543,7 @@ function ODAEProblem{iip}(sys,
     u0 = ModelingToolkit.varmap_to_vars(u0map, dvs; defaults = defs, tofloat = true)
     p = ModelingToolkit.varmap_to_vars(parammap, ps; defaults = defs, tofloat, use_union)
 
-    has_difference = any(isdifferenceeq, eqs)
-    cbs = process_events(sys; callback, has_difference, kwargs...)
+    cbs = process_events(sys; callback, kwargs...)
 
     kwargs = filter_kwargs(kwargs)
     if cbs === nothing
diff --git a/src/systems/callbacks.jl b/src/systems/callbacks.jl
index b97112e9e8..98bfc6a86a 100644
--- a/src/systems/callbacks.jl
+++ b/src/systems/callbacks.jl
@@ -498,7 +498,7 @@ merge_cb(::Nothing, x) = merge_cb(x, nothing)
 merge_cb(x, ::Nothing) = x
 merge_cb(x, y) = CallbackSet(x, y)
 
-function process_events(sys; callback = nothing, has_difference = false, kwargs...)
+function process_events(sys; callback = nothing, kwargs...)
     if has_continuous_events(sys)
         contin_cb = generate_rootfinding_callback(sys; kwargs...)
     else
@@ -509,9 +509,6 @@ function process_events(sys; callback = nothing, has_difference = false, kwargs.
     else
         discrete_cb = nothing
     end
-    difference_cb = has_difference ? generate_difference_cb(sys; kwargs...) : nothing
 
-    cb = merge_cb(contin_cb, difference_cb)
-    cb = merge_cb(cb, callback)
-    (discrete_cb === nothing) ? cb : CallbackSet(cb, discrete_cb...)
+    (discrete_cb === nothing) ? contin_cb : CallbackSet(contin_cb, discrete_cb...)
 end
diff --git a/src/systems/diffeqs/abstractodesystem.jl b/src/systems/diffeqs/abstractodesystem.jl
index 88a1ed9bb3..6df6d6fb09 100644
--- a/src/systems/diffeqs/abstractodesystem.jl
+++ b/src/systems/diffeqs/abstractodesystem.jl
@@ -144,7 +144,6 @@ function generate_function(sys::AbstractODESystem, dvs = states(sys), ps = param
         ddvs = implicit_dae ? map(Differential(get_iv(sys)), dvs) :
                nothing,
         isdde = false,
-        has_difference = false,
         kwargs...)
     if isdde
         eqs = delay_to_function(sys)
@@ -167,8 +166,7 @@ function generate_function(sys::AbstractODESystem, dvs = states(sys), ps = param
     if isdde
         build_function(rhss, u, DDE_HISTORY_FUN, p, t; kwargs...)
     else
-        pre, sol_states = get_substitutions_and_solved_states(sys,
-            no_postprocess = has_difference)
+        pre, sol_states = get_substitutions_and_solved_states(sys)
 
         if implicit_dae
             build_function(rhss, ddvs, u, p, t; postprocess_fbody = pre,
@@ -226,52 +224,6 @@ function delay_to_function(expr, iv, sts, ps, h)
     end
 end
 
-function generate_difference_cb(sys::ODESystem, dvs = states(sys), ps = parameters(sys);
-        kwargs...)
-    eqs = equations(sys)
-    check_operator_variables(eqs, Difference)
-
-    var2eq = Dict(arguments(eq.lhs)[1] => eq for eq in eqs if isdifference(eq.lhs))
-
-    u = map(x -> time_varying_as_func(value(x), sys), dvs)
-    p = map(x -> time_varying_as_func(value(x), sys), ps)
-    t = get_iv(sys)
-
-    body = map(dvs) do v
-        eq = get(var2eq, v, nothing)
-        eq === nothing && return v
-        d = operation(eq.lhs)
-        d.update ? eq.rhs : eq.rhs + v
-    end
-
-    pre = get_postprocess_fbody(sys)
-    cpre = get_preprocess_constants(body)
-    pre2 = x -> pre(cpre(x))
-    f_oop, f_iip = build_function(body, u, p, t; expression = Val{false},
-        postprocess_fbody = pre2, kwargs...)
-
-    cb_affect! = let f_oop = f_oop, f_iip = f_iip
-        function cb_affect!(integ)
-            if DiffEqBase.isinplace(integ.sol.prob)
-                tmp, = DiffEqBase.get_tmp_cache(integ)
-                f_iip(tmp, integ.u, integ.p, integ.t) # aliasing `integ.u` would be bad.
-                copyto!(integ.u, tmp)
-            else
-                integ.u = f_oop(integ.u, integ.p, integ.t)
-            end
-            return nothing
-        end
-    end
-
-    getdt(eq) = operation(eq.lhs).dt
-    deqs = values(var2eq)
-    dt = getdt(first(deqs))
-    all(dt == getdt(eq) for eq in deqs) ||
-        error("All difference variables should have same time steps.")
-
-    PeriodicCallback(cb_affect!, first(dt))
-end
-
 function calculate_massmatrix(sys::AbstractODESystem; simplify = false)
     eqs = [eq for eq in equations(sys) if !isdifferenceeq(eq)]
     dvs = states(sys)
@@ -940,11 +892,11 @@ function DiffEqBase.ODEProblem{iip, specialize}(sys::AbstractODESystem, u0map =
         check_length = true,
         kwargs...) where {iip, specialize}
     has_difference = any(isdifferenceeq, equations(sys))
+    has_difference && error("The operators Difference and DiscreteUpdate are deprecated. Use ShiftIndex instead.")
     f, u0, p = process_DEProblem(ODEFunction{iip, specialize}, sys, u0map, parammap;
         t = tspan !== nothing ? tspan[1] : tspan,
-        has_difference = has_difference,
         check_length, kwargs...)
-    cbs = process_events(sys; callback, has_difference, kwargs...)
+    cbs = process_events(sys; callback, kwargs...)
     if has_discrete_subsystems(sys) && (dss = get_discrete_subsystems(sys)) !== nothing
         affects, clocks, svs = ModelingToolkit.generate_discrete_affect(dss...)
         discrete_cbs = map(affects, clocks, svs) do affect, clock, sv
@@ -1003,23 +955,17 @@ function DiffEqBase.DAEProblem{iip}(sys::AbstractODESystem, du0map, u0map, tspan
         parammap = DiffEqBase.NullParameters();
         check_length = true, kwargs...) where {iip}
     has_difference = any(isdifferenceeq, equations(sys))
+    has_difference && error("The operators Difference and DiscreteUpdate are deprecated. Use ShiftIndex instead.")
     f, du0, u0, p = process_DEProblem(DAEFunction{iip}, sys, u0map, parammap;
-        implicit_dae = true, du0map = du0map,
-        has_difference = has_difference, check_length,
+        implicit_dae = true, du0map = du0map, check_length,
         kwargs...)
     diffvars = collect_differential_variables(sys)
     sts = states(sys)
     differential_vars = map(Base.Fix2(in, diffvars), sts)
     kwargs = filter_kwargs(kwargs)
 
-    if has_difference
-        DAEProblem{iip}(f, du0, u0, tspan, p;
-            difference_cb = generate_difference_cb(sys; kwargs...),
-            differential_vars = differential_vars, kwargs...)
-    else
-        DAEProblem{iip}(f, du0, u0, tspan, p; differential_vars = differential_vars,
+    DAEProblem{iip}(f, du0, u0, tspan, p; differential_vars = differential_vars,
             kwargs...)
-    end
 end
 
 function generate_history(sys::AbstractODESystem, u0; kwargs...)
@@ -1036,15 +982,15 @@ function DiffEqBase.DDEProblem{iip}(sys::AbstractODESystem, u0map = [],
         check_length = true,
         kwargs...) where {iip}
     has_difference = any(isdifferenceeq, equations(sys))
+    has_difference && error("The operators Difference and DiscreteUpdate are deprecated. Use ShiftIndex instead.")
     f, u0, p = process_DEProblem(DDEFunction{iip}, sys, u0map, parammap;
         t = tspan !== nothing ? tspan[1] : tspan,
-        has_difference = has_difference,
         symbolic_u0 = true,
         check_length, kwargs...)
     h_oop, h_iip = generate_history(sys, u0)
     h = h_oop
     u0 = h(p, tspan[1])
-    cbs = process_events(sys; callback, has_difference, kwargs...)
+    cbs = process_events(sys; callback, kwargs...)
     if has_discrete_subsystems(sys) && (dss = get_discrete_subsystems(sys)) !== nothing
         affects, clocks, svs = ModelingToolkit.generate_discrete_affect(dss...)
         discrete_cbs = map(affects, clocks, svs) do affect, clock, sv
@@ -1089,16 +1035,16 @@ function DiffEqBase.SDDEProblem{iip}(sys::AbstractODESystem, u0map = [],
         sparsenoise = nothing,
         kwargs...) where {iip}
     has_difference = any(isdifferenceeq, equations(sys))
+    has_difference && error("The operators Difference and DiscreteUpdate are deprecated. Use ShiftIndex instead.")
     f, u0, p = process_DEProblem(SDDEFunction{iip}, sys, u0map, parammap;
         t = tspan !== nothing ? tspan[1] : tspan,
-        has_difference = has_difference,
         symbolic_u0 = true,
         check_length, kwargs...)
     h_oop, h_iip = generate_history(sys, u0)
     h(out, p, t) = h_iip(out, p, t)
     h(p, t) = h_oop(p, t)
     u0 = h(p, tspan[1])
-    cbs = process_events(sys; callback, has_difference, kwargs...)
+    cbs = process_events(sys; callback, kwargs...)
     if has_discrete_subsystems(sys) && (dss = get_discrete_subsystems(sys)) !== nothing
         affects, clocks, svs = ModelingToolkit.generate_discrete_affect(dss...)
         discrete_cbs = map(affects, clocks, svs) do affect, clock, sv
diff --git a/src/systems/discrete_system/discrete_system.jl b/src/systems/discrete_system/discrete_system.jl
deleted file mode 100644
index dedb7bf6b3..0000000000
--- a/src/systems/discrete_system/discrete_system.jl
+++ /dev/null
@@ -1,464 +0,0 @@
-"""
-$(TYPEDEF)
-
-A system of difference equations.
-
-# Fields
-$(FIELDS)
-
-# Example
-
-```
-using ModelingToolkit
-
-@parameters σ=28.0 ρ=10.0 β=8/3 δt=0.1
-@variables t x(t)=1.0 y(t)=0.0 z(t)=0.0
-D = Difference(t; dt=δt)
-
-eqs = [D(x) ~ σ*(y-x),
-       D(y) ~ x*(ρ-z)-y,
-       D(z) ~ x*y - β*z]
-
-@named de = DiscreteSystem(eqs,t,[x,y,z],[σ,ρ,β]; tspan = (0, 1000.0)) # or
-@named de = DiscreteSystem(eqs)
-```
-"""
-struct DiscreteSystem <: AbstractTimeDependentSystem
-    """
-    A tag for the system. If two systems have the same tag, then they are
-    structurally identical.
-    """
-    tag::UInt
-    """The differential equations defining the discrete system."""
-    eqs::Vector{Equation}
-    """Independent variable."""
-    iv::BasicSymbolic{Real}
-    """Dependent (state) variables. Must not contain the independent variable."""
-    states::Vector
-    """Parameter variables. Must not contain the independent variable."""
-    ps::Vector
-    """Time span."""
-    tspan::Union{NTuple{2, Any}, Nothing}
-    """Array variables."""
-    var_to_name::Any
-    """Control parameters (some subset of `ps`)."""
-    ctrls::Vector
-    """Observed states."""
-    observed::Vector{Equation}
-    """
-    The name of the system
-    """
-    name::Symbol
-    """
-    The internal systems. These are required to have unique names.
-    """
-    systems::Vector{DiscreteSystem}
-    """
-    The default values to use when initial conditions and/or
-    parameters are not supplied in `DiscreteProblem`.
-    """
-    defaults::Dict
-    """
-    Inject assignment statements before the evaluation of the RHS function.
-    """
-    preface::Any
-    """
-    Type of the system.
-    """
-    connector_type::Any
-    """
-    Metadata for the system, to be used by downstream packages.
-    """
-    metadata::Any
-    """
-    Metadata for MTK GUI.
-    """
-    gui_metadata::Union{Nothing, GUIMetadata}
-    """
-    Cache for intermediate tearing state.
-    """
-    tearing_state::Any
-    """
-    Substitutions generated by tearing.
-    """
-    substitutions::Any
-    """
-    If a model `sys` is complete, then `sys.x` no longer performs namespacing.
-    """
-    complete::Bool
-    """
-    The hierarchical parent system before simplification.
-    """
-    parent::Any
-
-    function DiscreteSystem(tag, discreteEqs, iv, dvs, ps, tspan, var_to_name, ctrls,
-            observed,
-            name,
-            systems, defaults, preface, connector_type,
-            metadata = nothing, gui_metadata = nothing,
-            tearing_state = nothing, substitutions = nothing,
-            complete = false, parent = nothing; checks::Union{Bool, Int} = true)
-        if checks == true || (checks & CheckComponents) > 0
-            check_variables(dvs, iv)
-            check_parameters(ps, iv)
-        end
-        if checks == true || (checks & CheckUnits) > 0
-            u = __get_unit_type(dvs, ps, iv, ctrls)
-            check_units(u, discreteEqs)
-        end
-        new(tag, discreteEqs, iv, dvs, ps, tspan, var_to_name, ctrls, observed, name,
-            systems,
-            defaults,
-            preface, connector_type, metadata, gui_metadata,
-            tearing_state, substitutions, complete, parent)
-    end
-end
-
-"""
-    $(TYPEDSIGNATURES)
-
-Constructs a DiscreteSystem.
-"""
-function DiscreteSystem(eqs::AbstractVector{<:Equation}, iv, dvs, ps;
-        controls = Num[],
-        observed = Num[],
-        systems = DiscreteSystem[],
-        tspan = nothing,
-        name = nothing,
-        default_u0 = Dict(),
-        default_p = Dict(),
-        defaults = _merge(Dict(default_u0), Dict(default_p)),
-        preface = nothing,
-        connector_type = nothing,
-        metadata = nothing,
-        gui_metadata = nothing,
-        kwargs...)
-    name === nothing &&
-        throw(ArgumentError("The `name` keyword must be provided. Please consider using the `@named` macro"))
-    eqs = scalarize(eqs)
-    iv′ = value(iv)
-    dvs′ = value.(dvs)
-    ps′ = value.(ps)
-    ctrl′ = value.(controls)
-
-    if !(isempty(default_u0) && isempty(default_p))
-        Base.depwarn("`default_u0` and `default_p` are deprecated. Use `defaults` instead.",
-            :DiscreteSystem, force = true)
-    end
-    defaults = todict(defaults)
-    defaults = Dict(value(k) => value(v) for (k, v) in pairs(defaults))
-
-    var_to_name = Dict()
-    process_variables!(var_to_name, defaults, dvs′)
-    process_variables!(var_to_name, defaults, ps′)
-    isempty(observed) || collect_var_to_name!(var_to_name, (eq.lhs for eq in observed))
-
-    sysnames = nameof.(systems)
-    if length(unique(sysnames)) != length(sysnames)
-        throw(ArgumentError("System names must be unique."))
-    end
-    DiscreteSystem(Threads.atomic_add!(SYSTEM_COUNT, UInt(1)),
-        eqs, iv′, dvs′, ps′, tspan, var_to_name, ctrl′, observed, name, systems,
-        defaults, preface, connector_type, metadata, gui_metadata, kwargs...)
-end
-
-function DiscreteSystem(eqs, iv = nothing; kwargs...)
-    eqs = scalarize(eqs)
-    # NOTE: this assumes that the order of algebraic equations doesn't matter
-    diffvars = OrderedSet()
-    allstates = OrderedSet()
-    ps = OrderedSet()
-    # reorder equations such that it is in the form of `diffeq, algeeq`
-    diffeq = Equation[]
-    algeeq = Equation[]
-    # initial loop for finding `iv`
-    if iv === nothing
-        for eq in eqs
-            if !(eq.lhs isa Number) # assume eq.lhs is either Differential or Number
-                iv = iv_from_nested_difference(eq.lhs)
-                break
-            end
-        end
-    end
-    iv = value(iv)
-    iv === nothing && throw(ArgumentError("Please pass in independent variables."))
-    for eq in eqs
-        collect_vars_difference!(allstates, ps, eq.lhs, iv)
-        collect_vars_difference!(allstates, ps, eq.rhs, iv)
-        if isdifferenceeq(eq)
-            diffvar, _ = var_from_nested_difference(eq.lhs)
-            isequal(iv, iv_from_nested_difference(eq.lhs)) ||
-                throw(ArgumentError("A DiscreteSystem can only have one independent variable."))
-            diffvar in diffvars &&
-                throw(ArgumentError("The difference variable $diffvar is not unique in the system of equations."))
-            push!(diffvars, diffvar)
-            push!(diffeq, eq)
-        else
-            push!(algeeq, eq)
-        end
-    end
-    algevars = setdiff(allstates, diffvars)
-    # the orders here are very important!
-    return DiscreteSystem(append!(diffeq, algeeq), iv,
-        collect(Iterators.flatten((diffvars, algevars))), ps; kwargs...)
-end
-
-"""
-    $(TYPEDSIGNATURES)
-
-Generates an DiscreteProblem from an DiscreteSystem.
-"""
-function SciMLBase.DiscreteProblem(sys::DiscreteSystem, u0map = [], tspan = get_tspan(sys),
-        parammap = SciMLBase.NullParameters();
-        eval_module = @__MODULE__,
-        eval_expression = true,
-        use_union = false,
-        kwargs...)
-    dvs = states(sys)
-    ps = parameters(sys)
-    eqs = equations(sys)
-    eqs = linearize_eqs(sys, eqs)
-    iv = get_iv(sys)
-
-    defs = defaults(sys)
-    defs = mergedefaults(defs, parammap, ps)
-    defs = mergedefaults(defs, u0map, dvs)
-
-    u0 = varmap_to_vars(u0map, dvs; defaults = defs, tofloat = false)
-    p = varmap_to_vars(parammap, ps; defaults = defs, tofloat = false, use_union)
-
-    rhss = [eq.rhs for eq in eqs]
-    u = dvs
-
-    f_gen = generate_function(sys; expression = Val{eval_expression},
-        expression_module = eval_module)
-    f_oop, _ = (drop_expr(@RuntimeGeneratedFunction(eval_module, ex)) for ex in f_gen)
-    f(u, p, iv) = f_oop(u, p, iv)
-    fd = DiscreteFunction(f; syms = Symbol.(dvs), indepsym = Symbol(iv),
-        paramsyms = Symbol.(ps), sys = sys)
-    DiscreteProblem(fd, u0, tspan, p; kwargs...)
-end
-
-function linearize_eqs(sys, eqs = get_eqs(sys); return_max_delay = false)
-    unique_states = unique(operation.(states(sys)))
-    max_delay = Dict(v => 0.0 for v in unique_states)
-
-    r = @rule ~t::(t -> istree(t) && any(isequal(operation(t)), operation.(states(sys))) && is_delay_var(get_iv(sys), t)) => begin
-        delay = get_delay_val(get_iv(sys), first(arguments(~t)))
-        if delay > max_delay[operation(~t)]
-            max_delay[operation(~t)] = delay
-        end
-        nothing
-    end
-    SymbolicUtils.Postwalk(r).(rhss(eqs))
-
-    if any(values(max_delay) .> 0)
-        dts = Dict(v => Any[] for v in unique_states)
-        state_ops = Dict(v => Any[] for v in unique_states)
-        for v in unique_states
-            for eq in eqs
-                if isdifferenceeq(eq) && istree(arguments(eq.lhs)[1]) &&
-                   isequal(v, operation(arguments(eq.lhs)[1]))
-                    append!(dts[v], [operation(eq.lhs).dt])
-                    append!(state_ops[v], [operation(eq.lhs)])
-                end
-            end
-        end
-
-        all(length.(unique.(values(state_ops))) .<= 1) ||
-            error("Each state should be used with single difference operator.")
-
-        dts_gcd = Dict()
-        for v in keys(dts)
-            dts_gcd[v] = (length(dts[v]) > 0) ? first(dts[v]) : nothing
-        end
-
-        lin_eqs = [v(get_iv(sys) - (t)) ~ v(get_iv(sys) - (t - dts_gcd[v]))
-                   for v in unique_states if max_delay[v] > 0 && dts_gcd[v] !== nothing
-                   for t in collect(max_delay[v]:(-dts_gcd[v]):0)[1:(end - 1)]]
-        eqs = vcat(eqs, lin_eqs)
-    end
-    if return_max_delay
-        return eqs, max_delay
-    end
-    eqs
-end
-
-function get_delay_val(iv, x)
-    delay = x - iv
-    isequal(delay > 0, true) && error("Forward delay not permitted")
-    return -delay
-end
-
-function generate_function(sys::DiscreteSystem, dvs = states(sys), ps = parameters(sys);
-        kwargs...)
-    eqs = equations(sys)
-    check_operator_variables(eqs, Difference)
-    rhss = [eq.rhs for eq in eqs]
-
-    u = map(x -> time_varying_as_func(value(x), sys), dvs)
-    p = map(x -> time_varying_as_func(value(x), sys), ps)
-    t = get_iv(sys)
-
-    build_function(rhss, u, p, t; kwargs...)
-    pre, sol_states = get_substitutions_and_solved_states(sys)
-    build_function(rhss, u, p, t; postprocess_fbody = pre, states = sol_states, kwargs...)
-end
-
-"""
-```julia
-SciMLBase.DiscreteFunction{iip}(sys::DiscreteSystem, dvs = states(sys),
-                                ps = parameters(sys);
-                                version = nothing,
-                                kwargs...) where {iip}
-```
-
-Create a `DiscreteFunction` from the [`DiscreteSystem`](@ref). The arguments
-`dvs` and `ps` are used to set the order of the dependent variable and parameter
-vectors, respectively.
-"""
-function SciMLBase.DiscreteFunction(sys::DiscreteSystem, args...; kwargs...)
-    DiscreteFunction{true}(sys, args...; kwargs...)
-end
-
-function SciMLBase.DiscreteFunction{true}(sys::DiscreteSystem, args...; kwargs...)
-    DiscreteFunction{true, SciMLBase.AutoSpecialize}(sys, args...; kwargs...)
-end
-
-function SciMLBase.DiscreteFunction{false}(sys::DiscreteSystem, args...; kwargs...)
-    DiscreteFunction{false, SciMLBase.FullSpecialize}(sys, args...; kwargs...)
-end
-
-function SciMLBase.DiscreteFunction{iip, specialize}(sys::DiscreteSystem,
-        dvs = states(sys),
-        ps = parameters(sys),
-        u0 = nothing;
-        version = nothing,
-        p = nothing,
-        t = nothing,
-        eval_expression = true,
-        eval_module = @__MODULE__,
-        analytic = nothing,
-        simplify = false,
-        kwargs...) where {iip, specialize}
-    f_gen = generate_function(sys, dvs, ps; expression = Val{eval_expression},
-        expression_module = eval_module, kwargs...)
-    f_oop, f_iip = eval_expression ?
-                   (drop_expr(@RuntimeGeneratedFunction(eval_module, ex)) for ex in f_gen) :
-                   f_gen
-    f(u, p, t) = f_oop(u, p, t)
-    f(du, u, p, t) = f_iip(du, u, p, t)
-
-    if specialize === SciMLBase.FunctionWrapperSpecialize && iip
-        if u0 === nothing || p === nothing || t === nothing
-            error("u0, p, and t must be specified for FunctionWrapperSpecialize on DiscreteFunction.")
-        end
-        f = SciMLBase.wrapfun_iip(f, (u0, u0, p, t))
-    end
-
-    observedfun = let sys = sys, dict = Dict()
-        function generate_observed(obsvar, u, p, t)
-            obs = get!(dict, value(obsvar)) do
-                build_explicit_observed_function(sys, obsvar)
-            end
-            obs(u, p, t)
-        end
-    end
-
-    DiscreteFunction{iip, specialize}(f;
-        sys = sys,
-        syms = Symbol.(states(sys)),
-        indepsym = Symbol(get_iv(sys)),
-        paramsyms = Symbol.(ps),
-        observed = observedfun,
-        analytic = analytic)
-end
-
-"""
-```julia
-DiscreteFunctionExpr{iip}(sys::DiscreteSystem, dvs = states(sys),
-                          ps = parameters(sys);
-                          version = nothing,
-                          kwargs...) where {iip}
-```
-
-Create a Julia expression for an `DiscreteFunction` from the [`DiscreteSystem`](@ref).
-The arguments `dvs` and `ps` are used to set the order of the dependent
-variable and parameter vectors, respectively.
-"""
-struct DiscreteFunctionExpr{iip} end
-struct DiscreteFunctionClosure{O, I} <: Function
-    f_oop::O
-    f_iip::I
-end
-(f::DiscreteFunctionClosure)(u, p, t) = f.f_oop(u, p, t)
-(f::DiscreteFunctionClosure)(du, u, p, t) = f.f_iip(du, u, p, t)
-
-function DiscreteFunctionExpr{iip}(sys::DiscreteSystem, dvs = states(sys),
-        ps = parameters(sys), u0 = nothing;
-        version = nothing, p = nothing,
-        linenumbers = false,
-        simplify = false,
-        kwargs...) where {iip}
-    f_oop, f_iip = generate_function(sys, dvs, ps; expression = Val{true}, kwargs...)
-
-    fsym = gensym(:f)
-    _f = :($fsym = $DiscreteFunctionClosure($f_oop, $f_iip))
-
-    ex = quote
-        $_f
-        DiscreteFunction{$iip}($fsym,
-            syms = $(Symbol.(states(sys))),
-            indepsym = $(QuoteNode(Symbol(get_iv(sys)))),
-            paramsyms = $(Symbol.(parameters(sys))))
-    end
-    !linenumbers ? striplines(ex) : ex
-end
-
-function DiscreteFunctionExpr(sys::DiscreteSystem, args...; kwargs...)
-    DiscreteFunctionExpr{true}(sys, args...; kwargs...)
-end
-
-function process_DiscreteProblem(constructor, sys::DiscreteSystem, u0map, parammap;
-        version = nothing,
-        linenumbers = true, parallel = SerialForm(),
-        eval_expression = true,
-        use_union = false,
-        tofloat = !use_union,
-        kwargs...)
-    eqs = equations(sys)
-    dvs = states(sys)
-    ps = parameters(sys)
-
-    u0, p, defs = get_u0_p(sys, u0map, parammap; tofloat, use_union)
-
-    check_eqs_u0(eqs, dvs, u0; kwargs...)
-
-    f = constructor(sys, dvs, ps, u0;
-        linenumbers = linenumbers, parallel = parallel,
-        syms = Symbol.(dvs), paramsyms = Symbol.(ps),
-        eval_expression = eval_expression, kwargs...)
-    return f, u0, p
-end
-
-function DiscreteProblemExpr(sys::DiscreteSystem, args...; kwargs...)
-    DiscreteProblemExpr{true}(sys, args...; kwargs...)
-end
-
-function DiscreteProblemExpr{iip}(sys::DiscreteSystem, u0map, tspan,
-        parammap = DiffEqBase.NullParameters();
-        check_length = true,
-        kwargs...) where {iip}
-    f, u0, p = process_DiscreteProblem(DiscreteFunctionExpr{iip}, sys, u0map, parammap;
-        check_length, kwargs...)
-    linenumbers = get(kwargs, :linenumbers, true)
-
-    ex = quote
-        f = $f
-        u0 = $u0
-        p = $p
-        tspan = $tspan
-        DiscreteProblem(f, u0, tspan, p; $(filter_kwargs(kwargs)...))
-    end
-    !linenumbers ? striplines(ex) : ex
-end
diff --git a/src/systems/validation.jl b/src/systems/validation.jl
index 90b0745cd2..af09953491 100644
--- a/src/systems/validation.jl
+++ b/src/systems/validation.jl
@@ -51,7 +51,6 @@ function get_unit(x::Union{Symbolics.ArrayOp, Symbolics.Arr, Symbolics.CallWithM
     get_literal_unit(x)
 end
 get_unit(op::Differential, args) = get_unit(args[1]) / get_unit(op.x)
-get_unit(op::Difference, args) = get_unit(args[1]) / get_unit(op.t)
 get_unit(op::typeof(getindex), args) = get_unit(args[1])
 get_unit(x::SciMLBase.NullParameters) = unitless
 get_unit(op::typeof(instream), args) = get_unit(args[1])
diff --git a/src/utils.jl b/src/utils.jl
index 19668689a1..b08f708339 100644
--- a/src/utils.jl
+++ b/src/utils.jl
@@ -1,5 +1,4 @@
 get_iv(D::Differential) = D.x
-get_iv(D::Difference) = D.t
 
 function make_operation(@nospecialize(op), args)
     if op === (*)
@@ -184,7 +183,7 @@ function check_equations(eqs, iv)
     end
 end
 """
-Get all the independent variables with respect to which differentials/differences are taken.
+Get all the independent variables with respect to which differentials are taken.
 """
 function collect_ivs_from_nested_operator!(ivs, x, target_op)
     if !istree(x)
@@ -195,10 +194,8 @@ function collect_ivs_from_nested_operator!(ivs, x, target_op)
         push!(ivs, get_iv(op))
         x = if target_op <: Differential
             op.x
-        elseif target_op <: Difference
-            op.t
         else
-            error("Unknown target op type in collect_ivs $target_op. Pass Difference or Differential")
+            error("Unknown target op type in collect_ivs $target_op. Pass Differential")
         end
         collect_ivs_from_nested_operator!(ivs, x, target_op)
     end
@@ -261,7 +258,7 @@ Throw error when difference/derivative operation occurs in the R.H.S.
 """
 @noinline function throw_invalid_operator(opvar, eq, op::Type)
     if op === Difference
-        optext = "difference"
+        error("The Difference operator is deprecated, use ShiftIndex instead")        
     elseif op === Differential
         optext = "derivative"
     end
@@ -293,8 +290,7 @@ function check_operator_variables(eqs, op::T) where {T}
         if length(tmp) == 1
             x = only(tmp)
             if op === Differential
-                # Having a difference is fine for ODEs
-                is_tmp_fine = isdifferential(x) || isdifference(x)
+                is_tmp_fine = isdifferential(x)
             else
                 is_tmp_fine = istree(x) && !(operation(x) isa op)
             end
@@ -319,23 +315,6 @@ isoperator(op) = expr -> isoperator(expr, op)
 isdifferential(expr) = isoperator(expr, Differential)
 isdiffeq(eq) = isdifferential(eq.lhs)
 
-isdifference(expr) = isoperator(expr, Difference)
-isdifferenceeq(eq) = isdifference(eq.lhs)
-
-function iv_from_nested_difference(x::Symbolic)
-    istree(x) || return x
-    operation(x) isa Difference ? iv_from_nested_difference(arguments(x)[1]) :
-    arguments(x)[1]
-end
-iv_from_nested_difference(x) = nothing
-
-var_from_nested_difference(x, i = 0) = (nothing, nothing)
-function var_from_nested_difference(x::Symbolic, i = 0)
-    istree(x) && operation(x) isa Difference ?
-    var_from_nested_difference(arguments(x)[1], i + 1) :
-    (x, i)
-end
-
 isvariable(x::Num)::Bool = isvariable(value(x))
 function isvariable(x)::Bool
     x isa Symbolic || return false
@@ -350,7 +329,6 @@ end
 Return a `Set` containing all variables in `x` that appear in
 
   - differential equations if `op = Differential`
-  - difference equations if `op = Differential`
 
 Example:
 
@@ -390,9 +368,6 @@ function vars!(vars, O; op = Differential)
     return vars
 end
 
-difference_vars(x) = vars(x; op = Difference)
-difference_vars!(vars, O) = vars!(vars, O; op = Difference)
-
 function collect_operator_variables(sys::AbstractSystem, args...)
     collect_operator_variables(equations(sys), args...)
 end
@@ -404,7 +379,7 @@ end
     collect_operator_variables(eqs::AbstractVector{Equation}, op)
 
 Return a `Set` containing all variables that have Operator `op` applied to them.
-See also [`collect_differential_variables`](@ref), [`collect_difference_variables`](@ref).
+See also [`collect_differential_variables`](@ref).
 """
 function collect_operator_variables(eqs::AbstractVector{Equation}, op)
     vars = Set()
@@ -420,7 +395,6 @@ function collect_operator_variables(eqs::AbstractVector{Equation}, op)
     return diffvars
 end
 collect_differential_variables(sys) = collect_operator_variables(sys, Differential)
-collect_difference_variables(sys) = collect_operator_variables(sys, Difference)
 
 """
     collect_applied_operators(x, op)
@@ -471,20 +445,6 @@ function collect_vars!(states, parameters, expr, iv)
     return nothing
 end
 
-function collect_vars_difference!(states, parameters, expr, iv)
-    if issym(expr)
-        collect_var!(states, parameters, expr, iv)
-    else
-        for var in vars(expr)
-            if istree(var) && operation(var) isa Difference
-                var, _ = var_from_nested_difference(var)
-            end
-            collect_var!(states, parameters, var, iv)
-        end
-    end
-    return nothing
-end
-
 function collect_var!(states, parameters, var, iv)
     isequal(var, iv) && return nothing
     if isparameter(var) || (istree(var) && isparameter(operation(var)))
diff --git a/test/discretesystem.jl b/test/discretesystem.jl
deleted file mode 100644
index e8b2a0114d..0000000000
--- a/test/discretesystem.jl
+++ /dev/null
@@ -1,206 +0,0 @@
-# Example: Compartmental models in epidemiology
-#=
-- https://github.com/epirecipes/sir-julia/blob/master/markdown/function_map/function_map.md
-- https://en.wikipedia.org/wiki/Compartmental_models_in_epidemiology#Deterministic_versus_stochastic_epidemic_models
-=#
-using ModelingToolkit, Test
-using ModelingToolkit: get_metadata
-
-@inline function rate_to_proportion(r, t)
-    1 - exp(-r * t)
-end;
-
-# Independent and dependent variables and parameters
-@parameters t c nsteps δt β γ
-@constants h = 1
-D = Difference(t; dt = 0.1)
-@variables S(t) I(t) R(t)
-
-infection = rate_to_proportion(β * c * I / (S * h + I + R), δt * h) * S
-recovery = rate_to_proportion(γ * h, δt) * I
-
-# Equations
-eqs = [D(S) ~ S - infection * h,
-    D(I) ~ I + infection - recovery,
-    D(R) ~ R + recovery]
-
-# System
-@named sys = DiscreteSystem(eqs, t, [S, I, R], [c, nsteps, δt, β, γ]; controls = [β, γ])
-syss = structural_simplify(sys)
-@test syss == syss
-
-for df in [
-    DiscreteFunction(sys, [S, I, R], [c, nsteps, δt, β, γ]),
-    eval(DiscreteFunctionExpr(sys, [S, I, R], [c, nsteps, δt, β, γ])),
-]
-
-    # iip
-    du = zeros(3)
-    u = collect(1:3)
-    p = collect(1:5)
-    df.f(du, u, p, 0)
-    @test du ≈ [0.01831563888873422, 0.9816849729159067, 4.999999388195359]
-
-    # oop
-    @test df.f(u, p, 0) ≈ [0.01831563888873422, 0.9816849729159067, 4.999999388195359]
-end
-
-# Problem
-u0 = [S => 990.0, I => 10.0, R => 0.0]
-p = [β => 0.05, c => 10.0, γ => 0.25, δt => 0.1, nsteps => 400]
-tspan = (0.0, ModelingToolkit.value(substitute(nsteps, p))) # value function (from Symbolics) is used to convert a Num to Float64
-prob_map = DiscreteProblem(sys, u0, tspan, p)
-@test prob_map.f.sys === sys
-
-# Solution
-using OrdinaryDiffEq
-sol_map = solve(prob_map, FunctionMap());
-@test sol_map[S] isa Vector
-
-# Using defaults constructor
-@parameters t c=10.0 nsteps=400 δt=0.1 β=0.05 γ=0.25
-Diff = Difference(t; dt = 0.1)
-@variables S(t)=990.0 I(t)=10.0 R(t)=0.0
-
-infection2 = rate_to_proportion(β * c * I / (S + I + R), δt) * S
-recovery2 = rate_to_proportion(γ, δt) * I
-
-eqs2 = [D(S) ~ S - infection2,
-    D(I) ~ I + infection2 - recovery2,
-    D(R) ~ R + recovery2]
-
-@named sys = DiscreteSystem(eqs2; controls = [β, γ], tspan)
-@test ModelingToolkit.defaults(sys) != Dict()
-
-prob_map2 = DiscreteProblem(sys)
-sol_map2 = solve(prob_map, FunctionMap());
-
-@test sol_map.u == sol_map2.u
-@test sol_map.prob.p == sol_map2.prob.p
-
-# Direct Implementation
-
-function sir_map!(u_diff, u, p, t)
-    (S, I, R) = u
-    (β, c, γ, δt) = p
-    N = S + I + R
-    infection = rate_to_proportion(β * c * I / N, δt) * S
-    recovery = rate_to_proportion(γ, δt) * I
-    @inbounds begin
-        u_diff[1] = S - infection
-        u_diff[2] = I + infection - recovery
-        u_diff[3] = R + recovery
-    end
-    nothing
-end;
-u0 = [990.0, 10.0, 0.0];
-p = [0.05, 10.0, 0.25, 0.1];
-prob_map = DiscreteProblem(sir_map!, u0, tspan, p);
-sol_map2 = solve(prob_map, FunctionMap());
-
-@test Array(sol_map) ≈ Array(sol_map2)
-
-# Delayed difference equation
-@parameters t
-@variables x(..) y(..) z(t)
-D1 = Difference(t; dt = 1.5)
-D2 = Difference(t; dt = 2)
-
-@test ModelingToolkit.is_delay_var(Symbolics.value(t), Symbolics.value(x(t - 2)))
-@test ModelingToolkit.is_delay_var(Symbolics.value(t), Symbolics.value(y(t - 1)))
-@test !ModelingToolkit.is_delay_var(Symbolics.value(t), Symbolics.value(z))
-@test_throws ErrorException ModelingToolkit.get_delay_val(Symbolics.value(t),
-    Symbolics.arguments(Symbolics.value(x(t +
-                                          2)))[1])
-@test_throws ErrorException z(t)
-
-# Equations
-eqs = [
-    D1(x(t)) ~ 0.4x(t) + 0.3x(t - 1.5) + 0.1x(t - 3),
-    D2(y(t)) ~ 0.3y(t) + 0.7y(t - 2) + 0.1z * h,
-]
-
-# System
-@named sys = DiscreteSystem(eqs, t, [x(t), x(t - 1.5), x(t - 3), y(t), y(t - 2), z], [])
-
-eqs2, max_delay = ModelingToolkit.linearize_eqs(sys; return_max_delay = true)
-
-@test max_delay[Symbolics.operation(Symbolics.value(x(t)))] ≈ 3
-@test max_delay[Symbolics.operation(Symbolics.value(y(t)))] ≈ 2
-
-linearized_eqs = [eqs
-    x(t - 3.0) ~ x(t - 1.5)
-    x(t - 1.5) ~ x(t)
-    y(t - 2.0) ~ y(t)]
-@test all(eqs2 .== linearized_eqs)
-
-# observed variable handling
-@variables t x(t) RHS(t)
-@parameters τ
-@named fol = DiscreteSystem([D(x) ~ (1 - x) / τ]; observed = [RHS ~ (1 - x) / τ * h])
-@test isequal(RHS, @nonamespace fol.RHS)
-RHS2 = RHS
-@unpack RHS = fol
-@test isequal(RHS, RHS2)
-
-@testset "Preface tests" begin
-    @parameters t
-    using OrdinaryDiffEq
-    using Symbolics
-    using DiffEqBase: isinplace
-    using ModelingToolkit
-    using SymbolicUtils.Code
-    using SymbolicUtils: Sym
-
-    c = [0]
-    f = function f(c, d::Vector{Float64}, u::Vector{Float64}, p, t::Float64, dt::Float64)
-        c .= [c[1] + 1]
-        d .= randn(length(u))
-        nothing
-    end
-
-    dummy_identity(x, _) = x
-    @register_symbolic dummy_identity(x, y)
-
-    u0 = ones(5)
-    p0 = Float64[]
-    syms = [Symbol(:a, i) for i in 1:5]
-    syms_p = Symbol[]
-    dt = 0.1
-    @assert isinplace(f, 6)
-    wf = let c = c, buffer = similar(u0), u = similar(u0), p = similar(p0), dt = dt
-        t -> (f(c, buffer, u, p, t, dt); buffer)
-    end
-
-    num = hash(f) ⊻ length(u0) ⊻ length(p0)
-    buffername = Symbol(:fmi_buffer_, num)
-
-    Δ = DiscreteUpdate(t; dt = dt)
-    us = map(s -> (@variables $s(t))[1], syms)
-    ps = map(s -> (@variables $s(t))[1], syms_p)
-    buffer, = @variables $buffername[1:length(u0)]
-    dummy_var = Sym{Any}(:_) # this is safe because _ cannot be a rvalue in Julia
-
-    ss = Iterators.flatten((us, ps))
-    vv = Iterators.flatten((u0, p0))
-    defs = Dict{Any, Any}(s => v for (s, v) in zip(ss, vv))
-
-    preface = [Assignment(dummy_var, SetArray(true, term(getfield, wf, Meta.quot(:u)), us))
-        Assignment(dummy_var, SetArray(true, term(getfield, wf, Meta.quot(:p)), ps))
-        Assignment(buffer, term(wf, t))]
-    eqs = map(1:length(us)) do i
-        Δ(us[i]) ~ dummy_identity(buffer[i], us[i])
-    end
-
-    @named sys = DiscreteSystem(eqs, t, us, ps; defaults = defs, preface = preface)
-    prob = DiscreteProblem(sys, [], (0.0, 1.0))
-    sol = solve(prob, FunctionMap(); dt = dt)
-    @test c[1] + 1 == length(sol)
-end
-
-@parameters t
-@variables x(t) y(t)
-D = Difference(t; dt = 0.1)
-testdict = Dict([:test => 1])
-@named sys = DiscreteSystem([D(x) ~ 1.0]; metadata = testdict)
-@test get_metadata(sys) == testdict
diff --git a/test/odesystem.jl b/test/odesystem.jl
index 9bcbbfc130..5bbfa5e046 100644
--- a/test/odesystem.jl
+++ b/test/odesystem.jl
@@ -501,50 +501,6 @@ sol = solve(prob, Tsit5())
     map((x, y) -> x[2] .+ 2y, sol[x], sol[y]),
     map((x, y) -> x[3] .+ 3y, sol[x], sol[y]))
 
-# Mixed Difference Differential equations
-@parameters t a b c d
-@variables x(t) y(t)
-δ = Differential(t)
-Δ = Difference(t; dt = 0.1)
-U = DiscreteUpdate(t; dt = 0.1)
-eqs = [δ(x) ~ a * x - b * x * y
-    δ(y) ~ -c * y + d * x * y
-    Δ(x) ~ y
-    U(y) ~ x + 1]
-@named de = ODESystem(eqs, t, [x, y], [a, b, c, d])
-@test generate_difference_cb(de) isa ModelingToolkit.DiffEqCallbacks.DiscreteCallback
-
-# doesn't work with ODEFunction
-# prob = ODEProblem(ODEFunction{false}(de),[1.0,1.0],(0.0,1.0),[1.5,1.0,3.0,1.0])
-
-prob = ODEProblem(de, [1.0, 1.0], (0.0, 1.0), [1.5, 1.0, 3.0, 1.0], check_length = false)
-@test prob.kwargs[:callback] isa ModelingToolkit.DiffEqCallbacks.DiscreteCallback
-
-sol = solve(prob, Tsit5(); callback = prob.kwargs[:callback],
-    tstops = prob.tspan[1]:0.1:prob.tspan[2], verbose = false)
-
-# Direct implementation
-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
-
-prob2 = ODEProblem(lotka, [1.0, 1.0], (0.0, 1.0), [1.5, 1.0, 3.0, 1.0])
-function periodic_difference_affect!(int)
-    int.u = [int.u[1] + int.u[2], int.u[1] + 1]
-    return nothing
-end
-
-difference_cb = ModelingToolkit.PeriodicCallback(periodic_difference_affect!, 0.1)
-
-sol2 = solve(prob2, Tsit5(); callback = difference_cb,
-    tstops = collect(prob.tspan[1]:0.1:prob.tspan[2])[2:end], verbose = false)
-
-@test_broken sol(0:0.01:1)[x] ≈ sol2(0:0.01:1)[1, :]
-@test_broken sol(0:0.01:1)[y] ≈ sol2(0:0.01:1)[2, :]
-
 using ModelingToolkit
 
 function submodel(; name)
diff --git a/test/runtests.jl b/test/runtests.jl
index e9b392dd4a..730aac766f 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -21,7 +21,6 @@ end
         @safetestset "Linearization Tests" include("linearize.jl")
         @safetestset "Input Output Test" include("input_output_handling.jl")
         @safetestset "Clock Test" include("clock.jl")
-        @safetestset "DiscreteSystem Test" include("discretesystem.jl")
         @safetestset "ODESystem Test" include("odesystem.jl")
         @safetestset "Dynamic Quantities Test" include("dq_units.jl")
         @safetestset "Unitful Quantities Test" include("units.jl")
diff --git a/test/units.jl b/test/units.jl
index c518af8459..d1a6c57640 100644
--- a/test/units.jl
+++ b/test/units.jl
@@ -100,16 +100,6 @@ D = Differential(t)
 eqs = D.(x) .~ v
 ODESystem(eqs, name = :sys)
 
-# Difference equation
-@parameters t [unit = u"s"] a [unit = u"s"^-1]
-@variables x(t) [unit = u"kg"]
-δ = Differential(t)
-D = Difference(t; dt = 0.1u"s")
-eqs = [
-    δ(x) ~ a * x,
-]
-de = ODESystem(eqs, t, [x], [a], name = :sys)
-
 # Nonlinear system
 @parameters a [unit = u"kg"^-1]
 @variables x [unit = u"kg"]
diff --git a/test/variable_utils.jl b/test/variable_utils.jl
index 3de486e6d8..85e079a68a 100644
--- a/test/variable_utils.jl
+++ b/test/variable_utils.jl
@@ -16,8 +16,7 @@ new = (((1 / β - 1) + δ) / γ)^(1 / (γ - 1))
 @test iszero(sol - new)
 
 # Continuous
-using ModelingToolkit: isdifferential, isdifference, vars, difference_vars,
-    collect_difference_variables, collect_differential_variables,
+using ModelingToolkit: isdifferential, vars, collect_differential_variables,
     collect_ivs
 @variables t u(t) y(t)
 D = Differential(t)
@@ -33,32 +32,3 @@ aov = ModelingToolkit.collect_applied_operators(eq, Differential)
 
 ts = collect_ivs([eq])
 @test ts == Set([t])
-
-# Discrete
-z = Difference(t; dt = 1, update = true)
-eq = z(y) ~ u
-v = difference_vars(eq)
-@test v == Set([z(y), u])
-
-ov = collect_difference_variables(eq)
-@test ov == Set(Any[y])
-
-aov = ModelingToolkit.collect_applied_operators(eq, Difference)
-@test aov == Set(Any[z(y)])
-
-ts = collect_ivs([eq], Difference)
-@test ts == Set([t])
-
-@variables t
-function UnitDelay(dt; name)
-    @variables u(t)=0.0 [input = true] y(t)=0.0 [output = true]
-    z = Difference(t; dt = dt, update = true)
-    eqs = [
-        z(y) ~ u,
-    ]
-    DiscreteSystem(eqs, t, name = name)
-end
-
-dt = 0.1
-@named int = UnitDelay(dt)
-collect_difference_variables(int) == Set(Any[@nonamespace(int.y)])