diff --git a/config/default_configs/default_config.yml b/config/default_configs/default_config.yml
index 20ec70dc9b5..bef48175fc8 100644
--- a/config/default_configs/default_config.yml
+++ b/config/default_configs/default_config.yml
@@ -66,7 +66,7 @@ max_newton_iters_ode:
   help: "Maximum number of Newton's method iterations (only for ODE algorithms that use Newton's method)"
   value: 1
 ode_algo:
-  help: "ODE algorithm [`ARS343` (default), `SSP333`, `IMKG343a`, `ODE.Euler`, `ODE.IMEXEuler`, `ODE.Rosenbrock23`, etc.]"
+  help: "ODE algorithm [`ARS343` (default), `SSP333`, `IMKG343a`, etc.]"
   value: "ARS343"
 krylov_rtol:
   help: "Relative tolerance of the Krylov method (only for ClimaTimeSteppers.jl; only used if `use_krylov_method` is `true`)"
diff --git a/src/callbacks/callbacks.jl b/src/callbacks/callbacks.jl
index ff13a84e647..e865faa2549 100644
--- a/src/callbacks/callbacks.jl
+++ b/src/callbacks/callbacks.jl
@@ -64,7 +64,7 @@ NVTX.@annotate function turb_conv_affect_filter!(integrator)
         TC.affect_filter!(edmf, grid, state, tc_params, t)
     end
 
-    # We're lying to OrdinaryDiffEq.jl, in order to avoid
+    # We're lying to SciMLBase.jl, in order to avoid
     # paying for an additional `∑tendencies!` call, which is required
     # to support supplying a continuous representation of the
     # solution.
diff --git a/src/dycore_equations_deprecated/sgs_flux_tendencies.jl b/src/dycore_equations_deprecated/sgs_flux_tendencies.jl
index 9822e681c2b..2836f4e4c36 100644
--- a/src/dycore_equations_deprecated/sgs_flux_tendencies.jl
+++ b/src/dycore_equations_deprecated/sgs_flux_tendencies.jl
@@ -1,5 +1,4 @@
 using LinearAlgebra
-import OrdinaryDiffEq as ODE
 import Logging
 import TerminalLoggers
 import LinearAlgebra as LA
diff --git a/src/initial_conditions/InitialConditions.jl b/src/initial_conditions/InitialConditions.jl
index 29cef59fd7c..58191b02439 100644
--- a/src/initial_conditions/InitialConditions.jl
+++ b/src/initial_conditions/InitialConditions.jl
@@ -29,7 +29,6 @@ import ..Parameters as CAP
 import ..TurbulenceConvection as TC
 import Thermodynamics as TD
 import AtmosphericProfilesLibrary as APL
-import OrdinaryDiffEq as ODE
 import SciMLBase
 import Dierckx
 
diff --git a/src/parameterized_tendencies/radiation/radiation.jl b/src/parameterized_tendencies/radiation/radiation.jl
index b11cb245dc6..89ee4c1b459 100644
--- a/src/parameterized_tendencies/radiation/radiation.jl
+++ b/src/parameterized_tendencies/radiation/radiation.jl
@@ -4,7 +4,6 @@
 
 import ClimaComms
 import ClimaCore: Device, DataLayouts, Geometry, Spaces, Fields, Operators
-import OrdinaryDiffEq as ODE
 import Insolation
 import Thermodynamics as TD
 import .Parameters as CAP
diff --git a/src/prognostic_equations/implicit/schur_complement_W.jl b/src/prognostic_equations/implicit/schur_complement_W.jl
index 9ae95ed5bc3..7587b7acbd8 100644
--- a/src/prognostic_equations/implicit/schur_complement_W.jl
+++ b/src/prognostic_equations/implicit/schur_complement_W.jl
@@ -226,7 +226,7 @@ Note: The matrix S = A𝕄ρ Aρ𝕄 + A𝕄𝔼 A𝔼𝕄 + A𝕄𝕄 - I is th
 the large -I block in A.
 =#
 
-# Function required by OrdinaryDiffEq.jl
+# Function required by ClimaTimeSteppers.jl
 linsolve!(::Type{Val{:init}}, f, u0; kwargs...) = _linsolve!
 _linsolve!(x, A, b, update_matrix = false; kwargs...) =
     LinearAlgebra.ldiv!(x, A, b)
diff --git a/src/solver/solve.jl b/src/solver/solve.jl
index ea08f573eed..67c4f320431 100644
--- a/src/solver/solve.jl
+++ b/src/solver/solve.jl
@@ -1,5 +1,4 @@
 import ClimaTimeSteppers as CTS
-import OrdinaryDiffEq as ODE
 
 struct EfficiencyStats{TS <: Tuple, WT}
     tspan::TS
diff --git a/src/solver/type_getters.jl b/src/solver/type_getters.jl
index 89eb5c16183..855ef295aea 100644
--- a/src/solver/type_getters.jl
+++ b/src/solver/type_getters.jl
@@ -8,7 +8,6 @@ import ClimaAtmos.RRTMGPInterface as RRTMGPI
 import ClimaAtmos as CA
 import LinearAlgebra
 import ClimaCore.Fields
-import OrdinaryDiffEq as ODE
 import ClimaTimeSteppers as CTS
 import DiffEqCallbacks as DECB
 
@@ -340,29 +339,13 @@ is_explicit_CTS_algo_type(alg_or_tableau) =
 is_imex_CTS_algo_type(alg_or_tableau) =
     alg_or_tableau <: CTS.IMEXARKAlgorithmName
 
-is_implicit_type(::typeof(ODE.IMEXEuler)) = true
-is_implicit_type(alg_or_tableau) =
-    alg_or_tableau <: Union{
-        ODE.OrdinaryDiffEqImplicitAlgorithm,
-        ODE.OrdinaryDiffEqAdaptiveImplicitAlgorithm,
-    } || is_imex_CTS_algo_type(alg_or_tableau)
-
-is_ordinary_diffeq_newton(::typeof(ODE.IMEXEuler)) = true
-is_ordinary_diffeq_newton(alg_or_tableau) =
-    alg_or_tableau <: Union{
-        ODE.OrdinaryDiffEqNewtonAlgorithm,
-        ODE.OrdinaryDiffEqNewtonAdaptiveAlgorithm,
-    }
+is_implicit_type(alg_or_tableau) = is_imex_CTS_algo_type(alg_or_tableau)
 
 is_imex_CTS_algo(::CTS.IMEXAlgorithm) = true
 is_imex_CTS_algo(::SciMLBase.AbstractODEAlgorithm) = false
 
-is_implicit(::ODE.OrdinaryDiffEqImplicitAlgorithm) = true
-is_implicit(::ODE.OrdinaryDiffEqAdaptiveImplicitAlgorithm) = true
 is_implicit(ode_algo) = is_imex_CTS_algo(ode_algo)
 
-is_rosenbrock(::ODE.Rosenbrock23) = true
-is_rosenbrock(::ODE.Rosenbrock32) = true
 is_rosenbrock(::SciMLBase.AbstractODEAlgorithm) = false
 use_transform(ode_algo) =
     !(is_imex_CTS_algo(ode_algo) || is_rosenbrock(ode_algo))
@@ -409,28 +392,13 @@ Returns the ode algorithm
 =#
 function ode_configuration(::Type{FT}, parsed_args) where {FT}
     ode_name = parsed_args["ode_algo"]
-    alg_or_tableau = if startswith(ode_name, "ODE.")
-        @warn "apply_limiter flag is ignored for OrdinaryDiffEq algorithms"
-        getproperty(ODE, Symbol(split(ode_name, ".")[2]))
-    else
-        getproperty(CTS, Symbol(ode_name))
-    end
+    alg_or_tableau = getproperty(CTS, Symbol(ode_name))
     @info "Using ODE config: `$alg_or_tableau`"
 
     if is_explicit_CTS_algo_type(alg_or_tableau)
         return CTS.ExplicitAlgorithm(alg_or_tableau())
     elseif !is_implicit_type(alg_or_tableau)
         return alg_or_tableau()
-    elseif is_ordinary_diffeq_newton(alg_or_tableau)
-        if parsed_args["max_newton_iters_ode"] == 1
-            error("OridinaryDiffEq requires at least 2 Newton iterations")
-        end
-        # κ like a relative tolerance; its default value in ODE is 0.01
-        nlsolve = ODE.NLNewton(;
-            κ = parsed_args["max_newton_iters_ode"] == 2 ? Inf : 0.01,
-            max_iter = parsed_args["max_newton_iters_ode"],
-        )
-        return alg_or_tableau(; linsolve = linsolve!, nlsolve)
     elseif is_imex_CTS_algo_type(alg_or_tableau)
         newtons_method = CTS.NewtonsMethod(;
             max_iters = parsed_args["max_newton_iters_ode"],
@@ -477,9 +445,6 @@ function get_callbacks(parsed_args, simulation, atmos, params)
     (; dt) = simulation
 
     callbacks = ()
-    if startswith(parsed_args["ode_algo"], "ODE.")
-        callbacks = (callbacks..., call_every_n_steps(dss_callback!))
-    end
     dt_save_to_disk = time_to_seconds(parsed_args["dt_save_to_disk"])
     if !(dt_save_to_disk == Inf)
         callbacks = (