diff --git a/.buildkite/Manifest.toml b/.buildkite/Manifest.toml
index ed74a907a2..4b1b7aed9f 100644
--- a/.buildkite/Manifest.toml
+++ b/.buildkite/Manifest.toml
@@ -2,7 +2,7 @@
 
 julia_version = "1.10.4"
 manifest_format = "2.0"
-project_hash = "652f90825513bed89484e0ca25f9a53f61dedb8c"
+project_hash = "a3b8e1623c974758b823be61fb12f149202fadbb"
 
 [[deps.ADTypes]]
 git-tree-sha1 = "99a6f5d0ce1c7c6afdb759daa30226f71c54f6b0"
@@ -1744,6 +1744,28 @@ git-tree-sha1 = "a8b2d333cd90562b58b977b4033739360b37fb1f"
 uuid = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
 version = "6.87.0"
 
+[[deps.OrdinaryDiffEqCore]]
+deps = ["ADTypes", "Adapt", "ArrayInterface", "DataStructures", "DiffEqBase", "DocStringExtensions", "EnumX", "FastBroadcast", "FastClosures", "FillArrays", "FunctionWrappersWrappers", "InteractiveUtils", "LinearAlgebra", "Logging", "MacroTools", "MuladdMacro", "Polyester", "PrecompileTools", "Preferences", "RecursiveArrayTools", "Reexport", "SciMLBase", "SciMLOperators", "SciMLStructures", "SimpleUnPack", "Static", "StaticArrayInterface", "StaticArraysCore", "TruncatedStacktraces"]
+git-tree-sha1 = "bc70a38874c801549d415c16a23ec632ad91afd1"
+uuid = "bbf590c4-e513-4bbe-9b18-05decba2e5d8"
+version = "1.4.0"
+weakdeps = ["EnzymeCore"]
+
+    [deps.OrdinaryDiffEqCore.extensions]
+    OrdinaryDiffEqCoreEnzymeCoreExt = "EnzymeCore"
+
+[[deps.OrdinaryDiffEqSSPRK]]
+deps = ["DiffEqBase", "FastBroadcast", "MuladdMacro", "OrdinaryDiffEqCore", "Polyester", "PrecompileTools", "Preferences", "RecursiveArrayTools", "Reexport", "Static", "StaticArrays"]
+git-tree-sha1 = "7dbe4ac56f930df5e9abd003cedb54e25cbbea86"
+uuid = "669c94d9-1f4b-4b64-b377-1aa079aa2388"
+version = "1.2.0"
+
+[[deps.OrdinaryDiffEqTsit5]]
+deps = ["DiffEqBase", "FastBroadcast", "LinearAlgebra", "MuladdMacro", "OrdinaryDiffEqCore", "PrecompileTools", "Preferences", "RecursiveArrayTools", "Reexport", "Static", "TruncatedStacktraces"]
+git-tree-sha1 = "96552f7d4619fabab4038a29ed37dd55e9eb513a"
+uuid = "b1df2697-797e-41e3-8120-5422d3b24e4a"
+version = "1.1.0"
+
 [[deps.P11Kit_jll]]
 deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
 git-tree-sha1 = "2cd396108e178f3ae8dedbd8e938a18726ab2fbf"
diff --git a/.buildkite/Project.toml b/.buildkite/Project.toml
index 99c47af1a9..281d085e2d 100644
--- a/.buildkite/Project.toml
+++ b/.buildkite/Project.toml
@@ -35,6 +35,8 @@ NCDatasets = "85f8d34a-cbdd-5861-8df4-14fed0d494ab"
 NVTX = "5da4648a-3479-48b8-97b9-01cb529c0a1f"
 OrderedCollections = "bac558e1-5e72-5ebc-8fee-abe8a469f55d"
 OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
+OrdinaryDiffEqSSPRK = "669c94d9-1f4b-4b64-b377-1aa079aa2388"
+OrdinaryDiffEqTsit5 = "b1df2697-797e-41e3-8120-5422d3b24e4a"
 PProf = "e4faabce-9ead-11e9-39d9-4379958e3056"
 Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
 PrettyTables = "08abe8d2-0d0c-5749-adfa-8a2ac140af0d"
diff --git a/examples/bickleyjet/bickleyjet_cg.jl b/examples/bickleyjet/bickleyjet_cg.jl
index d772403d1e..574c250a69 100644
--- a/examples/bickleyjet/bickleyjet_cg.jl
+++ b/examples/bickleyjet/bickleyjet_cg.jl
@@ -12,7 +12,7 @@ import ClimaCore:
     Quadratures
 import ClimaCore.Geometry: ⊗
 
-using OrdinaryDiffEq: ODEProblem, solve, SSPRK33
+using OrdinaryDiffEqSSPRK: ODEProblem, solve, SSPRK33
 
 import Logging
 import TerminalLoggers
diff --git a/examples/bickleyjet/bickleyjet_cg_invariant_hypervisc.jl b/examples/bickleyjet/bickleyjet_cg_invariant_hypervisc.jl
index 6585c237b0..22f1be9cea 100644
--- a/examples/bickleyjet/bickleyjet_cg_invariant_hypervisc.jl
+++ b/examples/bickleyjet/bickleyjet_cg_invariant_hypervisc.jl
@@ -11,7 +11,7 @@ import ClimaCore:
     Topologies,
     Quadratures,
     DataLayouts
-using OrdinaryDiffEq: ODEProblem, solve, SSPRK33
+using OrdinaryDiffEqSSPRK: ODEProblem, solve, SSPRK33
 
 using Logging
 ClimaComms.@import_required_backends
diff --git a/examples/bickleyjet/bickleyjet_cg_unsmesh.jl b/examples/bickleyjet/bickleyjet_cg_unsmesh.jl
index 0d29c80b8e..7c33deab11 100644
--- a/examples/bickleyjet/bickleyjet_cg_unsmesh.jl
+++ b/examples/bickleyjet/bickleyjet_cg_unsmesh.jl
@@ -12,7 +12,7 @@ import ClimaCore:
     Quadratures
 import ClimaCore.Geometry: ⊗
 
-using OrdinaryDiffEq: ODEProblem, solve, SSPRK33
+using OrdinaryDiffEqSSPRK: ODEProblem, solve, SSPRK33
 
 import Logging
 import TerminalLoggers
diff --git a/examples/bickleyjet/bickleyjet_dg.jl b/examples/bickleyjet/bickleyjet_dg.jl
index 3e7978ff30..cf8ff15afa 100644
--- a/examples/bickleyjet/bickleyjet_dg.jl
+++ b/examples/bickleyjet/bickleyjet_dg.jl
@@ -14,7 +14,7 @@ import ClimaCore:
 import ClimaCore.Geometry: ⊗
 import ClimaCore.RecursiveApply: ⊞, rdiv, rmap
 
-using OrdinaryDiffEq: ODEProblem, solve, SSPRK33
+using OrdinaryDiffEqSSPRK: ODEProblem, solve, SSPRK33
 
 import Logging
 import TerminalLoggers
diff --git a/examples/column/advect.jl b/examples/column/advect.jl
index 7853ab976d..1ac06cd3c8 100644
--- a/examples/column/advect.jl
+++ b/examples/column/advect.jl
@@ -10,7 +10,7 @@ import ClimaCore:
     Geometry,
     Spaces
 
-using OrdinaryDiffEq: ODEProblem, solve, SSPRK33
+using OrdinaryDiffEqSSPRK: ODEProblem, solve, SSPRK33
 
 import Logging
 import TerminalLoggers
diff --git a/examples/column/advect_diffusion.jl b/examples/column/advect_diffusion.jl
index 5bee5e4467..4caa1444a3 100644
--- a/examples/column/advect_diffusion.jl
+++ b/examples/column/advect_diffusion.jl
@@ -10,7 +10,7 @@ import ClimaCore:
     Geometry,
     Spaces
 
-using OrdinaryDiffEq: ODEProblem, solve, SSPRK33
+using OrdinaryDiffEqSSPRK: ODEProblem, solve, SSPRK33
 
 import Logging
 import TerminalLoggers
diff --git a/examples/column/bb_fct_advection.jl b/examples/column/bb_fct_advection.jl
index 8c049679f6..cc0a823efb 100644
--- a/examples/column/bb_fct_advection.jl
+++ b/examples/column/bb_fct_advection.jl
@@ -2,7 +2,7 @@ import ClimaComms
 ClimaComms.@import_required_backends
 using Test
 using LinearAlgebra
-using OrdinaryDiffEq: ODEProblem, solve
+using SciMLBase: ODEProblem, solve
 using ClimaTimeSteppers
 
 import ClimaCore:
diff --git a/examples/column/ekman.jl b/examples/column/ekman.jl
index 37860cd63d..362b14a5ba 100644
--- a/examples/column/ekman.jl
+++ b/examples/column/ekman.jl
@@ -11,7 +11,7 @@ import ClimaCore:
     Geometry,
     Spaces
 
-using OrdinaryDiffEq: OrdinaryDiffEq, ODEProblem, solve, SSPRK33
+using OrdinaryDiffEqSSPRK: ODEProblem, solve, SSPRK33
 
 import Logging
 import TerminalLoggers
diff --git a/examples/column/fct_advection.jl b/examples/column/fct_advection.jl
index 6134ce4d41..fee5881916 100644
--- a/examples/column/fct_advection.jl
+++ b/examples/column/fct_advection.jl
@@ -2,7 +2,7 @@ import ClimaComms
 ClimaComms.@import_required_backends
 using Test
 using LinearAlgebra
-using OrdinaryDiffEq: ODEProblem, solve, SSPRK33
+using OrdinaryDiffEqSSPRK: ODEProblem, solve, SSPRK33
 
 import ClimaCore:
     Fields,
diff --git a/examples/column/heat.jl b/examples/column/heat.jl
index cabe614ac2..ddfda96e95 100644
--- a/examples/column/heat.jl
+++ b/examples/column/heat.jl
@@ -10,7 +10,7 @@ import ClimaCore:
     Geometry,
     Spaces
 
-using OrdinaryDiffEq: ODEProblem, solve, SSPRK33
+using OrdinaryDiffEqSSPRK: ODEProblem, solve, SSPRK33
 
 import Logging
 import TerminalLoggers
diff --git a/examples/column/hydrostatic.jl b/examples/column/hydrostatic.jl
index c2dd942b2e..e5aa39a354 100644
--- a/examples/column/hydrostatic.jl
+++ b/examples/column/hydrostatic.jl
@@ -10,7 +10,7 @@ import ClimaCore:
     Geometry,
     Spaces
 
-using OrdinaryDiffEq: OrdinaryDiffEq, ODEProblem, solve, SSPRK33
+using OrdinaryDiffEqSSPRK: ODEProblem, solve, SSPRK33
 
 import Logging
 import TerminalLoggers
diff --git a/examples/column/hydrostatic_discrete.jl b/examples/column/hydrostatic_discrete.jl
index f543ae6d46..5d3a2848f6 100644
--- a/examples/column/hydrostatic_discrete.jl
+++ b/examples/column/hydrostatic_discrete.jl
@@ -10,7 +10,7 @@ import ClimaCore:
     Geometry,
     Spaces
 
-using OrdinaryDiffEq: OrdinaryDiffEq, ODEProblem, solve, SSPRK33
+using OrdinaryDiffEqSSPRK: ODEProblem, solve, SSPRK33
 
 import Logging
 import TerminalLoggers
diff --git a/examples/column/hydrostatic_ekman.jl b/examples/column/hydrostatic_ekman.jl
index 483f037dd8..40aa393590 100644
--- a/examples/column/hydrostatic_ekman.jl
+++ b/examples/column/hydrostatic_ekman.jl
@@ -14,7 +14,7 @@ import ClimaCore:
 
 import ClimaCore.Geometry: ⊗
 
-using OrdinaryDiffEq: OrdinaryDiffEq, ODEProblem, solve, SSPRK33
+using OrdinaryDiffEqSSPRK: ODEProblem, solve, SSPRK33
 
 import Logging
 import TerminalLoggers
diff --git a/examples/column/step.jl b/examples/column/step.jl
index 70755a3eaa..08fbe186f1 100644
--- a/examples/column/step.jl
+++ b/examples/column/step.jl
@@ -11,7 +11,7 @@ import ClimaCore:
     Geometry,
     Spaces
 
-using OrdinaryDiffEq: ODEProblem, solve, SSPRK33
+using OrdinaryDiffEqSSPRK: ODEProblem, solve, SSPRK33
 
 import Logging
 import TerminalLoggers
diff --git a/examples/column/wave.jl b/examples/column/wave.jl
index 9573e1dba1..2cde0c40cd 100644
--- a/examples/column/wave.jl
+++ b/examples/column/wave.jl
@@ -11,7 +11,7 @@ import ClimaCore:
     Geometry,
     Spaces
 
-using OrdinaryDiffEq: OrdinaryDiffEq, ODEProblem, solve, SSPRK33
+using OrdinaryDiffEqSSPRK: ODEProblem, solve, SSPRK33
 
 import Logging
 import TerminalLoggers
diff --git a/examples/column/zalesak_fct_advection.jl b/examples/column/zalesak_fct_advection.jl
index 9dc13a86a7..a2e1dff9ab 100644
--- a/examples/column/zalesak_fct_advection.jl
+++ b/examples/column/zalesak_fct_advection.jl
@@ -2,7 +2,7 @@ import ClimaComms
 ClimaComms.@import_required_backends
 using Test
 using LinearAlgebra
-using OrdinaryDiffEq: ODEProblem, solve
+using SciMLBase: ODEProblem, solve
 using ClimaTimeSteppers
 
 import ClimaCore:
diff --git a/examples/hybrid/box/bubble_3d_invariant_rhoe.jl b/examples/hybrid/box/bubble_3d_invariant_rhoe.jl
index a6bbca5b87..4f2112e4e0 100644
--- a/examples/hybrid/box/bubble_3d_invariant_rhoe.jl
+++ b/examples/hybrid/box/bubble_3d_invariant_rhoe.jl
@@ -13,7 +13,8 @@ using Adapt
 using ClimaComms
 ClimaComms.@import_required_backends
 FloatType = eval(Meta.parse(get(ARGS, 1, "Float64")))
-using StaticArrays, IntervalSets, LinearAlgebra, OrdinaryDiffEq
+using StaticArrays, IntervalSets, LinearAlgebra, SciMLBase
+using OrdinaryDiffEqSSPRK: SSPRK33
 using DocStringExtensions
 
 import ClimaCore:
@@ -446,7 +447,7 @@ function bubble_3d_invariant_ρe(ARGS, comms_ctx, ::Type{FT}) where {FT}
     # run!
     Δt = sim_params.Δt
     prob = ODEProblem(rhs_invariant!, Y, (0.0, sim_params.t_int), ghost_buffer)
-    integrator = OrdinaryDiffEq.init(
+    integrator = SciMLBase.init(
         prob,
         SSPRK33(),
         dt = Δt,
@@ -460,7 +461,7 @@ function bubble_3d_invariant_ρe(ARGS, comms_ctx, ::Type{FT}) where {FT}
         throw(:exit_profile)
     end
 
-    t_diff = @elapsed sol_invariant = OrdinaryDiffEq.solve!(integrator)
+    t_diff = @elapsed sol_invariant = SciMLBase.solve!(integrator)
 
     if ClimaComms.iamroot(comms_ctx)
         println("Walltime = $t_diff seconds")
diff --git a/examples/hybrid/box/bubble_3d_invariant_rhotheta.jl b/examples/hybrid/box/bubble_3d_invariant_rhotheta.jl
index 1712ee8d77..4169b0838d 100644
--- a/examples/hybrid/box/bubble_3d_invariant_rhotheta.jl
+++ b/examples/hybrid/box/bubble_3d_invariant_rhotheta.jl
@@ -304,7 +304,7 @@ dYdt = similar(Y);
 rhs_invariant!(dYdt, Y, nothing, 0.0);
 
 # run!
-using OrdinaryDiffEq
+using OrdinaryDiffEqSSPRK: ODEProblem, solve, SSPRK33
 Δt = 0.050
 prob = ODEProblem(rhs_invariant!, Y, (0.0, 1.0))
 sol = solve(
diff --git a/examples/hybrid/box/bubble_3d_rhotheta.jl b/examples/hybrid/box/bubble_3d_rhotheta.jl
index ab32a7df03..0fa78f1ac1 100644
--- a/examples/hybrid/box/bubble_3d_rhotheta.jl
+++ b/examples/hybrid/box/bubble_3d_rhotheta.jl
@@ -305,10 +305,10 @@ rhs!(dYdt, Y, nothing, 0.0);
 
 
 # run!
-using OrdinaryDiffEq
+using OrdinaryDiffEqSSPRK: ODEProblem, init, solve!, SSPRK33
 Δt = 0.05
 prob = ODEProblem(rhs!, Y, (0.0, 1.0))
-integrator = OrdinaryDiffEq.init(
+integrator = init(
     prob,
     SSPRK33(),
     dt = Δt,
@@ -321,7 +321,7 @@ if haskey(ENV, "CI_PERF_SKIP_RUN") # for performance analysis
     throw(:exit_profile)
 end
 
-sol = @timev OrdinaryDiffEq.solve!(integrator)
+sol = @timev solve!(integrator)
 
 ENV["GKSwstype"] = "nul"
 import Plots
diff --git a/examples/hybrid/driver.jl b/examples/hybrid/driver.jl
index 6b72273018..c0347f408d 100644
--- a/examples/hybrid/driver.jl
+++ b/examples/hybrid/driver.jl
@@ -47,7 +47,7 @@ atexit() do
     global_logger(prev_logger)
 end
 
-using OrdinaryDiffEq
+using SciMLBase
 using DiffEqCallbacks
 using JLD2
 
@@ -163,7 +163,7 @@ problem = ODE.ODEProblem(
     (t_start, t_end),
     p,
 )
-integrator = OrdinaryDiffEq.init(
+integrator = SciMLBase.init(
     problem,
     ode_algo;
     saveat = dt_save_to_sol == 0 ? [] : dt_save_to_sol,
@@ -182,7 +182,7 @@ end
 @info "Running `$test_dir/$test_file_name` test case"
 @info "on a vertical $z_stretch_string grid"
 
-walltime = @elapsed sol = OrdinaryDiffEq.solve!(integrator)
+walltime = @elapsed sol = SciMLBase.solve!(integrator)
 any(isnan, sol.u[end]) && error("NaNs found in result.")
 
 if is_distributed # replace sol.u on the root processor with the global sol.u
diff --git a/examples/hybrid/plane/bubble_2d_invariant_rhoe.jl b/examples/hybrid/plane/bubble_2d_invariant_rhoe.jl
index 1842e4fe31..7c3c619186 100644
--- a/examples/hybrid/plane/bubble_2d_invariant_rhoe.jl
+++ b/examples/hybrid/plane/bubble_2d_invariant_rhoe.jl
@@ -253,10 +253,10 @@ dYdt = similar(Y);
 rhs_invariant!(dYdt, Y, nothing, 0.0);
 
 # run!
-using OrdinaryDiffEq
+using OrdinaryDiffEqSSPRK: init, ODEProblem, solve!, SSPRK33
 Δt = 0.04
 prob = ODEProblem(rhs_invariant!, Y, (0.0, 1200.0))
-integrator = OrdinaryDiffEq.init(
+integrator = init(
     prob,
     SSPRK33(),
     dt = Δt,
@@ -269,7 +269,7 @@ if haskey(ENV, "CI_PERF_SKIP_RUN") # for performance analysis
     throw(:exit_profile)
 end
 
-sol = @timev OrdinaryDiffEq.solve!(integrator)
+sol = @timev solve!(integrator)
 
 ENV["GKSwstype"] = "nul"
 import Plots, ClimaCorePlots
diff --git a/examples/hybrid/plane/density_current_2d_flux_form.jl b/examples/hybrid/plane/density_current_2d_flux_form.jl
index fed1d3257b..f2da579811 100644
--- a/examples/hybrid/plane/density_current_2d_flux_form.jl
+++ b/examples/hybrid/plane/density_current_2d_flux_form.jl
@@ -295,11 +295,11 @@ rhs!(dYdt, Y, nothing, 0.0);
 
 
 # run!
-using OrdinaryDiffEq
+using OrdinaryDiffEqSSPRK: ODEProblem, init, solve!, SSPRK33
 Δt = 0.3
 prob = ODEProblem(rhs!, Y, (0.0, 900.0))
 
-integrator = OrdinaryDiffEq.init(
+integrator = init(
     prob,
     SSPRK33(),
     dt = Δt,
@@ -312,7 +312,7 @@ if haskey(ENV, "CI_PERF_SKIP_RUN") # for performance analysis
     throw(:exit_profile)
 end
 
-sol = @timev OrdinaryDiffEq.solve!(integrator)
+sol = @timev solve!(integrator)
 
 ENV["GKSwstype"] = "nul"
 using ClimaCorePlots, Plots
diff --git a/examples/hybrid/plane/density_current_2dinvariant_rhoe.jl b/examples/hybrid/plane/density_current_2dinvariant_rhoe.jl
index 83b36c4066..43eb7a2862 100644
--- a/examples/hybrid/plane/density_current_2dinvariant_rhoe.jl
+++ b/examples/hybrid/plane/density_current_2dinvariant_rhoe.jl
@@ -272,11 +272,11 @@ dYdt = similar(Y);
 rhs_invariant!(dYdt, Y, nothing, 0.0);
 
 # run!
-using OrdinaryDiffEq
+using OrdinaryDiffEqSSPRK: ODEProblem, init, solve!, SSPRK33
 timeend = 900.0
 Δt = 0.3
 prob = ODEProblem(rhs_invariant!, Y, (0.0, timeend))
-integrator = OrdinaryDiffEq.init(
+integrator = init(
     prob,
     SSPRK33(),
     dt = Δt,
@@ -289,7 +289,7 @@ if haskey(ENV, "CI_PERF_SKIP_RUN") # for performance analysis
     throw(:exit_profile)
 end
 
-sol = @timev OrdinaryDiffEq.solve!(integrator)
+sol = @timev solve!(integrator)
 
 ENV["GKSwstype"] = "nul"
 import Plots, ClimaCorePlots
diff --git a/examples/hybrid/plane/isothermal_channel.jl b/examples/hybrid/plane/isothermal_channel.jl
index 65b1f99230..c700c3f3aa 100644
--- a/examples/hybrid/plane/isothermal_channel.jl
+++ b/examples/hybrid/plane/isothermal_channel.jl
@@ -277,10 +277,10 @@ dYdt = similar(Y);
 rhs_invariant!(dYdt, Y, nothing, 0.0);
 
 # run!
-using OrdinaryDiffEq
+using OrdinaryDiffEqSSPRK: ODEProblem, init, solve!, SSPRK33
 Δt = 0.5
 prob = ODEProblem(rhs_invariant!, Y, (0.0, 15000.0))
-integrator = OrdinaryDiffEq.init(
+integrator = init(
     prob,
     SSPRK33(),
     dt = Δt,
@@ -293,7 +293,7 @@ if haskey(ENV, "CI_PERF_SKIP_RUN") # for performance analysis
     throw(:exit_profile)
 end
 
-sol = @timev OrdinaryDiffEq.solve!(integrator)
+sol = @timev solve!(integrator)
 
 ENV["GKSwstype"] = "nul"
 import Plots, ClimaCorePlots
diff --git a/examples/hybrid/plane/topo_agnesi_nh.jl b/examples/hybrid/plane/topo_agnesi_nh.jl
index 2480bb320f..121fe79dfb 100644
--- a/examples/hybrid/plane/topo_agnesi_nh.jl
+++ b/examples/hybrid/plane/topo_agnesi_nh.jl
@@ -335,7 +335,7 @@ dYdt = similar(Y);
 rhs_invariant!(dYdt, Y, nothing, 0.0);
 
 # run!
-using OrdinaryDiffEq
+using OrdinaryDiffEqSSPRK: ODEProblem, init, solve!, SSPRK33
 Δt = 1.0
 timeend = 72000.0
 function make_dss_func()
@@ -347,7 +347,7 @@ end
 dss_func = make_dss_func()
 dss_callback = FunctionCallingCallback(dss_func, func_start = true)
 prob = ODEProblem(rhs_invariant!, Y, (0.0, timeend))
-integrator = OrdinaryDiffEq.init(
+integrator = init(
     prob,
     SSPRK33(),
     dt = Δt,
@@ -361,7 +361,7 @@ if haskey(ENV, "CI_PERF_SKIP_RUN") # for performance analysis
     throw(:exit_profile)
 end
 
-sol = @timev OrdinaryDiffEq.solve!(integrator)
+sol = @timev solve!(integrator)
 
 ENV["GKSwstype"] = "nul"
 import Plots, ClimaCorePlots
diff --git a/examples/hybrid/plane/topo_schar_nh.jl b/examples/hybrid/plane/topo_schar_nh.jl
index 9a5e905dec..a0f7324837 100644
--- a/examples/hybrid/plane/topo_schar_nh.jl
+++ b/examples/hybrid/plane/topo_schar_nh.jl
@@ -397,7 +397,7 @@ dYdt = similar(Y);
 rhs_invariant!(dYdt, Y, nothing, 0.0);
 
 # run!
-using OrdinaryDiffEq
+using OrdinaryDiffEqSSPRK: ODEProblem, init, solve!, SSPRK33
 Δt = min(Lx / nx / np / 300, Lz / nz / 300) * 0.50
 @show Δt
 
@@ -411,7 +411,7 @@ end
 dss_func = make_dss_func()
 dss_callback = FunctionCallingCallback(dss_func, func_start = true)
 prob = ODEProblem(rhs_invariant!, Y, (0.0, timeend))
-integrator = OrdinaryDiffEq.init(
+integrator = init(
     prob,
     SSPRK33(),
     dt = Δt,
@@ -425,7 +425,7 @@ if haskey(ENV, "CI_PERF_SKIP_RUN") # for performance analysis
     throw(:exit_profile)
 end
 
-sol = @timev OrdinaryDiffEq.solve!(integrator)
+sol = @timev solve!(integrator)
 
 ENV["GKSwstype"] = "nul"
 import Plots, ClimaCorePlots
diff --git a/examples/hybrid/sphere/deformation_flow.jl b/examples/hybrid/sphere/deformation_flow.jl
index 5e2999758d..9b4ce6e5c8 100644
--- a/examples/hybrid/sphere/deformation_flow.jl
+++ b/examples/hybrid/sphere/deformation_flow.jl
@@ -1,6 +1,6 @@
 import ClimaComms
 ClimaComms.@import_required_backends
-using OrdinaryDiffEq
+using SciMLBase: ODEProblem, init, solve
 using Test
 using Statistics: mean
 
diff --git a/examples/hybrid/sphere/hadley_circulation.jl b/examples/hybrid/sphere/hadley_circulation.jl
index b489a182b2..0ada1ad4cb 100644
--- a/examples/hybrid/sphere/hadley_circulation.jl
+++ b/examples/hybrid/sphere/hadley_circulation.jl
@@ -16,7 +16,7 @@ import ClimaCore:
     Operators
 import ClimaCore.Utilities: half
 
-using OrdinaryDiffEq: ODEProblem, solve
+using SciMLBase: ODEProblem, solve
 using ClimaTimeSteppers
 
 using NCDatasets, ClimaCoreTempestRemap
diff --git a/examples/hybrid/sphere/nonhydrostatic_gravity_wave.jl b/examples/hybrid/sphere/nonhydrostatic_gravity_wave.jl
index ba618afdf7..b5495bfb53 100644
--- a/examples/hybrid/sphere/nonhydrostatic_gravity_wave.jl
+++ b/examples/hybrid/sphere/nonhydrostatic_gravity_wave.jl
@@ -15,7 +15,7 @@ import ClimaCore:
     Fields,
     Operators
 
-using OrdinaryDiffEq: ODEProblem, solve, SSPRK33
+using OrdinaryDiffEqSSPRK: ODEProblem, solve, SSPRK33
 
 import Logging
 import TerminalLoggers
@@ -251,7 +251,7 @@ dYdt = similar(Y)
 rhs!(dYdt, Y, nothing, 0.0)
 
 # run!
-using OrdinaryDiffEq
+using OrdinaryDiffEqSSPRK: ODEProblem, solve, SSPRK33
 # Solve the ODE
 time_end = 600
 dt = 3
diff --git a/examples/hybrid/sphere/solid_body_rotation_3d.jl b/examples/hybrid/sphere/solid_body_rotation_3d.jl
index 2c9eae36ca..ee0fb5d191 100644
--- a/examples/hybrid/sphere/solid_body_rotation_3d.jl
+++ b/examples/hybrid/sphere/solid_body_rotation_3d.jl
@@ -15,7 +15,7 @@ import ClimaCore:
     Operators
 import ClimaCore.Utilities: half
 
-using OrdinaryDiffEq: ODEProblem, solve, SSPRK33
+using OrdinaryDiffEqSSPRK: ODEProblem, solve, SSPRK33
 
 import Logging
 import TerminalLoggers
@@ -279,13 +279,13 @@ dYdt = similar(Y)
 rhs!(dYdt, Y, parameters, 0.0)
 
 # run!
-using OrdinaryDiffEq
+using OrdinaryDiffEqSSPRK: ODEProblem, init, solve!, SSPRK33
 # Solve the ODE
 T = 3600
 dt = 5
 prob = ODEProblem(rhs!, Y, (0.0, T), parameters)
 
-integrator = OrdinaryDiffEq.init(
+integrator = init(
     prob,
     SSPRK33(),
     dt = dt,
@@ -300,7 +300,7 @@ if haskey(ENV, "CI_PERF_SKIP_RUN") # for performance analysis
 end
 
 # solve ode
-sol = @timev OrdinaryDiffEq.solve!(integrator)
+sol = @timev solve!(integrator)
 
 uₕ_phy = Geometry.transform.(Ref(Geometry.UVAxis()), sol.u[end].uₕ)
 w_phy = Geometry.transform.(Ref(Geometry.WAxis()), sol.u[end].w)
diff --git a/examples/ordinary_diff_eq_bug_fixes.jl b/examples/ordinary_diff_eq_bug_fixes.jl
index cbf22b9e68..0a818e93cd 100644
--- a/examples/ordinary_diff_eq_bug_fixes.jl
+++ b/examples/ordinary_diff_eq_bug_fixes.jl
@@ -13,10 +13,10 @@ using OrdinaryDiffEq:
     isJcurrent,
     postamble!,
     TryAgain
-using OrdinaryDiffEq.SciMLBase: RECOMPILE_BY_DEFAULT
+using SciMLBase: RECOMPILE_BY_DEFAULT
 
 import OrdinaryDiffEq: nlsolve!
-import OrdinaryDiffEq.SciMLBase: SplitFunction
+import SciMLBase: SplitFunction
 
 #=
 Issue: calc_W! attempts to access f.Wfact when has_Wfact(f) and f.Wfact_t when
diff --git a/examples/plane/limiters_advection.jl b/examples/plane/limiters_advection.jl
index 81d3b9f289..17f34ea927 100644
--- a/examples/plane/limiters_advection.jl
+++ b/examples/plane/limiters_advection.jl
@@ -12,7 +12,7 @@ import ClimaCore:
     Limiters,
     Quadratures
 
-using OrdinaryDiffEq: ODEProblem, solve
+using SciMLBase: ODEProblem, solve
 using ClimaTimeSteppers
 
 import Logging
diff --git a/examples/sphere/limiters_advection.jl b/examples/sphere/limiters_advection.jl
index 31295528ee..fd6b4e8e75 100644
--- a/examples/sphere/limiters_advection.jl
+++ b/examples/sphere/limiters_advection.jl
@@ -12,7 +12,8 @@ import ClimaCore:
     Limiters,
     Quadratures
 
-using OrdinaryDiffEq, Test
+using Test
+using OrdinaryDiffEqSSPRK: ODEProblem, solve, SSPRK33
 
 import Logging
 import TerminalLoggers
diff --git a/examples/sphere/shallow_water.jl b/examples/sphere/shallow_water.jl
index 8b3b3cb4d9..309e428859 100644
--- a/examples/sphere/shallow_water.jl
+++ b/examples/sphere/shallow_water.jl
@@ -18,8 +18,7 @@ import ClimaCore:
     DataLayouts
 
 import QuadGK
-import OrdinaryDiffEq
-using OrdinaryDiffEq: ODEProblem, solve, SSPRK33
+using OrdinaryDiffEqSSPRK: ODEProblem, init, solve!, SSPRK33
 
 using Logging
 using ClimaComms
@@ -582,7 +581,7 @@ function shallow_water_driver(ARGS, ::Type{FT}) where {FT}
     T = 60 * 60 * 24 * 2
 
     prob = ODEProblem(rhs!, Y, (0.0, T), parameters)
-    integrator = OrdinaryDiffEq.init(
+    integrator = init(
         prob,
         SSPRK33(),
         dt = dt,
@@ -593,10 +592,10 @@ function shallow_water_driver(ARGS, ::Type{FT}) where {FT}
     )
 
     if usempi
-        walltime = @elapsed sol = OrdinaryDiffEq.solve!(integrator)
+        walltime = @elapsed sol = solve!(integrator)
         ClimaComms.iamroot(context) && println("walltime = $walltime (sec)")
     else
-        sol = @timev OrdinaryDiffEq.solve!(integrator)
+        sol = @timev solve!(integrator)
     end
     sol_global = []
 
diff --git a/examples/sphere/shallow_water_cuda.jl b/examples/sphere/shallow_water_cuda.jl
index 36b3f6c259..a65a0d4fa7 100644
--- a/examples/sphere/shallow_water_cuda.jl
+++ b/examples/sphere/shallow_water_cuda.jl
@@ -4,7 +4,6 @@ using DocStringExtensions
 using LinearAlgebra
 using ClimaTimeSteppers
 import SciMLBase
-import OrdinaryDiffEq as ODE
 import ClimaTimeSteppers as CTS
 using DiffEqCallbacks
 using NVTX
@@ -568,13 +567,13 @@ function shallow_water_driver_cuda(ARGS, ::Type{FT}) where {FT}
         end
         return nothing
     end
-    prob = ODE.ODEProblem(
+    prob = SciMLBase.ODEProblem(
         CTS.ClimaODEFunction(; T_exp! = rhs!),
         Y,
         (FT(0), T),
         parameters,
     )
-    integrator = ODE.init(
+    integrator = SciMLBase.init(
         prob,
         CTS.ExplicitAlgorithm(CTS.SSP33ShuOsher()),
         dt = dt,
diff --git a/examples/sphere/solidbody.jl b/examples/sphere/solidbody.jl
index be1c03be98..73149d49ab 100644
--- a/examples/sphere/solidbody.jl
+++ b/examples/sphere/solidbody.jl
@@ -11,7 +11,7 @@ import ClimaCore:
     Topologies,
     Quadratures
 
-using OrdinaryDiffEq: ODEProblem, solve, SSPRK33
+using OrdinaryDiffEqSSPRK: ODEProblem, solve, SSPRK33
 
 import Logging
 import TerminalLoggers
diff --git a/perf/benchmark.jl b/perf/benchmark.jl
index 20185e3fcf..f891bfa45a 100644
--- a/perf/benchmark.jl
+++ b/perf/benchmark.jl
@@ -17,7 +17,8 @@ catch err
     end
 end
 
-OrdinaryDiffEq.step!(integrator) # compile first
-trial = BenchmarkTools.@benchmark OrdinaryDiffEq.step!($integrator)
+using SciMLBase: solve!, step!
+step!(integrator) # compile first
+trial = BenchmarkTools.@benchmark step!($integrator)
 show(stdout, MIME("text/plain"), trial)
 println()
diff --git a/perf/flame.jl b/perf/flame.jl
index 15fcda30d0..b45656b83d 100644
--- a/perf/flame.jl
+++ b/perf/flame.jl
@@ -17,11 +17,12 @@ catch err
     end
 end
 
-OrdinaryDiffEq.step!(integrator) # compile first
+using SciMLBase: step!, solve!
+step!(integrator) # compile first
 Profile.clear_malloc_data()
 prof = Profile.@profile begin
     for _ in 1:100
-        OrdinaryDiffEq.step!(integrator)
+        step!(integrator)
     end
 end
 
diff --git a/test/Operators/finitedifference/convergence_advection_diffusion1d.jl b/test/Operators/finitedifference/convergence_advection_diffusion1d.jl
index f67b665bc3..6dc561ec44 100644
--- a/test/Operators/finitedifference/convergence_advection_diffusion1d.jl
+++ b/test/Operators/finitedifference/convergence_advection_diffusion1d.jl
@@ -2,7 +2,7 @@ import ClimaComms
 ClimaComms.@import_required_backends
 using Test
 using LinearAlgebra
-using OrdinaryDiffEq: ODEProblem, solve, SSPRK33, Tsit5
+using OrdinaryDiffEqTsit5: ODEProblem, solve, Tsit5
 
 import ClimaCore:
     Fields,
diff --git a/test/Operators/hybrid/simulation_3d.jl b/test/Operators/hybrid/simulation_3d.jl
index 4244ba860f..1777fa1e4f 100644
--- a/test/Operators/hybrid/simulation_3d.jl
+++ b/test/Operators/hybrid/simulation_3d.jl
@@ -1,5 +1,6 @@
 include("utils_3d.jl")
 
+using OrdinaryDiffEqSSPRK: ODEProblem, solve, SSPRK33
 device = ClimaComms.device()
 
 @testset "2D SE, 1D FV Extruded Domain ∇ ODE Solve vertical" begin
@@ -25,7 +26,6 @@ device = ClimaComms.device()
     dudt = zeros(eltype(U), hv_center_space)
     rhs!(dudt, U, nothing, 0.0)
 
-    using OrdinaryDiffEq
     Δt = 0.01
     prob = ODEProblem(rhs!, U, (0.0, 2π))
     sol = solve(prob, SSPRK33(), dt = Δt)
@@ -60,7 +60,6 @@ end
     dudt = zeros(eltype(U), hv_center_space)
     rhs!(dudt, U, nothing, 0.0)
 
-    using OrdinaryDiffEq
     Δt = 0.01
     prob = ODEProblem(rhs!, U, (0.0, 2π))
     sol = solve(prob, SSPRK33(), dt = Δt)
@@ -140,7 +139,6 @@ end
         return h
     end
 
-    using OrdinaryDiffEq
     U = Fields.FieldVector(h = h_init(0.5, 0.5, 0.5))
     U_fieldbc = copy(U)
     Δt = 0.01
diff --git a/test/Operators/hybrid/simulation_cuda.jl b/test/Operators/hybrid/simulation_cuda.jl
index 60648515f1..b63534078b 100644
--- a/test/Operators/hybrid/simulation_cuda.jl
+++ b/test/Operators/hybrid/simulation_cuda.jl
@@ -1,5 +1,5 @@
 include("utils_cuda.jl")
-using OrdinaryDiffEq
+using OrdinaryDiffEqSSPRK: ODEProblem, solve, SSPRK33
 
 @testset "2D SE, 1D FD Extruded Domain ∇ ODE Solve horizontal CUDA" begin