diff --git a/src/ensemble/basic_ensemble_solve.jl b/src/ensemble/basic_ensemble_solve.jl
index e4dd54fac..c831679eb 100644
--- a/src/ensemble/basic_ensemble_solve.jl
+++ b/src/ensemble/basic_ensemble_solve.jl
@@ -105,7 +105,7 @@ function batch_func(i, prob, alg; kwargs...)
     new_prob = prob.prob_func(_prob, i, iter)
     rerun = true
     x = prob.output_func(solve(new_prob, alg; kwargs...), i)
-    if !(typeof(x) <: Tuple)
+    if !(x isa Tuple)
         rerun_warn()
         _x = (x, false)
     else
@@ -117,7 +117,7 @@ function batch_func(i, prob, alg; kwargs...)
         _prob = prob.safetycopy ? deepcopy(prob.prob) : prob.prob
         new_prob = prob.prob_func(_prob, i, iter)
         x = prob.output_func(solve(new_prob, alg; kwargs...), i)
-        if !(typeof(x) <: Tuple)
+        if !(x isa Tuple)
             rerun_warn()
             _x = (x, false)
         else
@@ -170,7 +170,7 @@ function solve_batch(prob, alg, ensemblealg::EnsembleThreads, II, pmap_batch_siz
         return solve_batch(prob, alg, EnsembleSerial(), II, pmap_batch_size; kwargs...)
     end
 
-    if typeof(prob.prob) <: AbstractJumpProblem && length(II) != 1
+    if prob.prob isa AbstractJumpProblem && length(II) != 1
         probs = [deepcopy(prob.prob) for i in 1:nthreads]
     else
         probs = prob.prob
diff --git a/src/ensemble/ensemble_analysis.jl b/src/ensemble/ensemble_analysis.jl
index 509869848..c68fd5dad 100644
--- a/src/ensemble/ensemble_analysis.jl
+++ b/src/ensemble/ensemble_analysis.jl
@@ -7,7 +7,7 @@ get_timestep(sim, i) = (getindex(sol, i) for sol in sim)
 get_timepoint(sim, t) = (sol(t) for sol in sim)
 function componentwise_vectors_timestep(sim, i)
     arr = [get_timestep(sim, i)...]
-    if typeof(arr[1]) <: AbstractArray
+    if arr[1] isa AbstractArray
         return vecarr_to_vectors(VectorOfArray(arr))
     else
         return arr
@@ -15,7 +15,7 @@ function componentwise_vectors_timestep(sim, i)
 end
 function componentwise_vectors_timepoint(sim, t)
     arr = [get_timepoint(sim, t)...]
-    if typeof(arr[1]) <: AbstractArray
+    if arr[1] isa AbstractArray
         return vecarr_to_vectors(VectorOfArray(arr))
     else
         return arr
@@ -123,7 +123,7 @@ end
 
 function SciMLBase.EnsembleSummary(sim::SciMLBase.AbstractEnsembleSolution{T, N},
     t = sim[1].t; quantiles = [0.05, 0.95]) where {T, N}
-    if typeof(sim[1]) <: SciMLSolution
+    if sim[1] isa SciMLSolution
         m, v = timeseries_point_meanvar(sim, t)
         med = timeseries_point_median(sim, t)
         qlow = timeseries_point_quantile(sim, quantiles[1], t)
@@ -190,13 +190,13 @@ function componentwise_mean(A)
     mean = zero(x0) ./ 1
     for x in A
         n += 1
-        if typeof(x0) <: AbstractArray && !(typeof(x0) <: StaticArraysCore.SArray)
+        if x0 isa AbstractArray && !(x0 isa StaticArraysCore.SArray)
             mean .+= x
         else
             mean += x
         end
     end
-    if typeof(x0) <: AbstractArray && !(typeof(x0) <: StaticArraysCore.SArray)
+    if x0 isa AbstractArray && !(x0 isa StaticArraysCore.SArray)
         mean ./= n
     else
         mean /= n
@@ -215,7 +215,7 @@ function componentwise_meanvar(A; bessel = true)
     delta2 = zero(x0) ./ 1
     for x in A
         n += 1
-        if typeof(x0) <: AbstractArray && !(typeof(x0) <: StaticArraysCore.SArray)
+        if x0 isa AbstractArray && !(x0 isa StaticArraysCore.SArray)
             delta .= x .- mean
             mean .+= delta ./ n
             delta2 .= x .- mean
@@ -231,13 +231,13 @@ function componentwise_meanvar(A; bessel = true)
         return NaN
     else
         if bessel
-            if typeof(x0) <: AbstractArray && !(typeof(x0) <: StaticArraysCore.SArray)
+            if x0 isa AbstractArray && !(x0 isa StaticArraysCore.SArray)
                 M2 .= M2 ./ (n .- 1)
             else
                 M2 = M2 ./ (n .- 1)
             end
         else
-            if typeof(x0) <: AbstractArray && !(typeof(x0) <: StaticArraysCore.SArray)
+            if x0 isa AbstractArray && !(x0 isa StaticArraysCore.SArray)
                 M2 .= M2 ./ n
             else
                 M2 = M2 ./ n
@@ -257,7 +257,7 @@ function componentwise_meancov(A, B; bessel = true)
     dx = zero(x0) ./ 1
     for (x, y) in zip(A, B)
         n += 1
-        if typeof(x0) <: AbstractArray && !(typeof(x0) <: StaticArraysCore.SArray)
+        if x0 isa AbstractArray && !(x0 isa StaticArraysCore.SArray)
             dx .= x .- meanx
             meanx .+= dx ./ n
             meany .+= (y .- meany) ./ n
@@ -273,13 +273,13 @@ function componentwise_meancov(A, B; bessel = true)
         return NaN
     else
         if bessel
-            if typeof(x0) <: AbstractArray && !(typeof(x0) <: StaticArraysCore.SArray)
+            if x0 isa AbstractArray && !(x0 isa StaticArraysCore.SArray)
                 C .= C ./ (n .- 1)
             else
                 C = C ./ (n .- 1)
             end
         else
-            if typeof(x0) <: AbstractArray && !(typeof(x0) <: StaticArraysCore.SArray)
+            if x0 isa AbstractArray && !(x0 isa StaticArraysCore.SArray)
                 C .= C ./ n
             else
                 C = C ./ n
@@ -293,7 +293,7 @@ function componentwise_meancor(A, B; bessel = true)
     mx, my, cov = componentwise_meancov(A, B; bessel = bessel)
     mx, vx = componentwise_meanvar(A; bessel = bessel)
     my, vy = componentwise_meanvar(B; bessel = bessel)
-    if typeof(vx) <: AbstractArray
+    if vx isa AbstractArray
         vx .= sqrt.(vx)
         vy .= sqrt.(vy)
     else
@@ -316,7 +316,7 @@ function componentwise_weighted_meancov(A, B, W; weight_type = :reliability)
     dx = zero(x0) ./ 1
     for (x, y, w) in zip(A, B, W)
         n += 1
-        if typeof(x0) <: AbstractArray && !(typeof(x0) <: StaticArraysCore.SArray)
+        if x0 isa AbstractArray && !(x0 isa StaticArraysCore.SArray)
             wsum .+= w
             wsum2 .+= w .* w
             dx .= x .- meanx
@@ -336,19 +336,19 @@ function componentwise_weighted_meancov(A, B, W; weight_type = :reliability)
         return NaN
     else
         if weight_type == :population
-            if typeof(x0) <: AbstractArray && !(typeof(x0) <: StaticArraysCore.SArray)
+            if x0 isa AbstractArray && !(x0 isa StaticArraysCore.SArray)
                 C .= C ./ wsum
             else
                 C = C ./ wsum
             end
         elseif weight_type == :reliability
-            if typeof(x0) <: AbstractArray && !(typeof(x0) <: StaticArraysCore.SArray)
+            if x0 isa AbstractArray && !(x0 isa StaticArraysCore.SArray)
                 C .= C ./ (wsum .- wsum2 ./ wsum)
             else
                 C = C ./ (wsum .- wsum2 ./ wsum)
             end
         elseif weight_type == :frequency
-            if typeof(x0) <: AbstractArray && !(typeof(x0) <: StaticArraysCore.SArray)
+            if x0 isa AbstractArray && !(x0 isa StaticArraysCore.SArray)
                 C .= C ./ (wsum .- 1)
             else
                 C = C ./ (wsum .- 1)
diff --git a/src/ensemble/ensemble_solutions.jl b/src/ensemble/ensemble_solutions.jl
index b0ad1d337..5bb198512 100644
--- a/src/ensemble/ensemble_solutions.jl
+++ b/src/ensemble/ensemble_solutions.jl
@@ -139,7 +139,7 @@ end
 ### Plot Recipes
 
 @recipe function f(sim::AbstractEnsembleSolution;
-    zcolors = typeof(sim.u) <: AbstractArray ? fill(nothing, length(sim.u)) :
+    zcolors = sim.u isa AbstractArray ? fill(nothing, length(sim.u)) :
               nothing,
     trajectories = eachindex(sim))
     for i in trajectories
@@ -156,16 +156,16 @@ end
 end
 
 @recipe function f(sim::EnsembleSummary;
-    trajectories = typeof(sim.u[1]) <: AbstractArray ? eachindex(sim.u[1]) :
+    trajectories = sim.u[1] isa AbstractArray ? eachindex(sim.u[1]) :
                    1,
     error_style = :ribbon, ci_type = :quantile)
     if ci_type == :SEM
-        if typeof(sim.u[1]) <: AbstractArray
+        if sim.u[1] isa AbstractArray
             u = vecarr_to_vectors(sim.u)
         else
             u = [sim.u.u]
         end
-        if typeof(sim.u[1]) <: AbstractArray
+        if sim.u[1] isa AbstractArray
             ci_low = vecarr_to_vectors(VectorOfArray([sqrt.(sim.v[i] / sim.num_monte) .*
                                                       1.96 for i in 1:length(sim.v)]))
             ci_high = ci_low
@@ -175,12 +175,12 @@ end
             ci_high = ci_low
         end
     elseif ci_type == :quantile
-        if typeof(sim.med[1]) <: AbstractArray
+        if sim.med[1] isa AbstractArray
             u = vecarr_to_vectors(sim.med)
         else
             u = [sim.med.u]
         end
-        if typeof(sim.u[1]) <: AbstractArray
+        if sim.u[1] isa AbstractArray
             ci_low = u - vecarr_to_vectors(sim.qlow)
             ci_high = vecarr_to_vectors(sim.qhigh) - u
         else
diff --git a/src/integrator_interface.jl b/src/integrator_interface.jl
index 52e6edb28..7046398a7 100644
--- a/src/integrator_interface.jl
+++ b/src/integrator_interface.jl
@@ -750,7 +750,7 @@ Base.length(iter::TimeChoiceIterator) = length(iter.ts)
 
 @recipe function f(integrator::DEIntegrator;
     denseplot = (integrator.opts.calck ||
-                 typeof(integrator) <: AbstractSDEIntegrator) &&
+                 integrator isa AbstractSDEIntegrator) &&
                 integrator.iter > 0,
     plotdensity = 10,
     plot_analytic = false, vars = nothing, idxs = nothing)
@@ -797,7 +797,7 @@ Base.length(iter::TimeChoiceIterator) = length(iter.ts)
             else # just get values
                 if x[j] == 0
                     push!(plot_vecs[j - 1], integrator.t)
-                elseif x[j] == 1 && !(typeof(integrator.u) <: AbstractArray)
+                elseif x[j] == 1 && !(integrator.u isa AbstractArray)
                     push!(plot_vecs[j - 1], integrator.u)
                 else
                     push!(plot_vecs[j - 1], integrator.u[x[j]])
@@ -816,7 +816,7 @@ Base.length(iter::TimeChoiceIterator) = length(iter.ts)
                 else # Just get values
                     if x[j] == 0
                         push!(plot_vecs[j], integrator.t)
-                    elseif x[j] == 1 && !(typeof(integrator.u) <: AbstractArray)
+                    elseif x[j] == 1 && !(integrator.u isa AbstractArray)
                         push!(plot_vecs[j],
                             integrator.sol.prob.f(Val{:analytic}, integrator.t,
                                 integrator.sol[1]))
@@ -840,7 +840,7 @@ Base.length(iter::TimeChoiceIterator) = length(iter.ts)
     end
 
     # Special case labels when idxs = (:x,:y,:z) or (:x) or [:x,:y] ...
-    if typeof(idxs) <: Tuple && (typeof(idxs[1]) == Symbol && typeof(idxs[2]) == Symbol)
+    if idxs isa Tuple && (typeof(idxs[1]) == Symbol && typeof(idxs[2]) == Symbol)
         xlabel --> idxs[1]
         ylabel --> idxs[2]
         if length(idxs) > 2
diff --git a/src/interpolation.jl b/src/interpolation.jl
index e64d91c35..d212ac82f 100644
--- a/src/interpolation.jl
+++ b/src/interpolation.jl
@@ -81,7 +81,7 @@ end
     continuity::Symbol = :left) where {I, D}
     t = id.t
     u = id.u
-    typeof(id) <: HermiteInterpolation && (du = id.du)
+    id isa HermiteInterpolation && (du = id.du)
     tdir = sign(t[end] - t[1])
     idx = sortperm(tvals, rev = tdir < 0)
     i = 2 # Start the search thinking it's between t[1] and t[2]
@@ -91,9 +91,9 @@ end
         error("Solution interpolation cannot extrapolate past the final timepoint. Either solve on a longer timespan or use the local extrapolation from the integrator interface.")
     tdir * tvals[idx[1]] < tdir * t[1] &&
         error("Solution interpolation cannot extrapolate before the first timepoint. Either start solving earlier or use the local extrapolation from the integrator interface.")
-    if typeof(idxs) <: Number
+    if idxs isa Number
         vals = Vector{eltype(first(u))}(undef, length(tvals))
-    elseif typeof(idxs) <: AbstractVector
+    elseif idxs isa AbstractVector
         vals = Vector{Vector{eltype(first(u))}}(undef, length(tvals))
     else
         vals = Vector{eltype(u)}(undef, length(tvals))
@@ -101,7 +101,7 @@ end
     for j in idx
         tval = tvals[j]
         i = searchsortedfirst(@view(t[i:end]), tval, rev = tdir < 0) + i - 1 # It's in the interval t[i-1] to t[i]
-        avoid_constant_ends = deriv != Val{0} #|| typeof(tval) <: ForwardDiff.Dual
+        avoid_constant_ends = deriv != Val{0} #|| tval isa ForwardDiff.Dual
         avoid_constant_ends && i == 1 && (i += 1)
         if !avoid_constant_ends && t[i - 1] == tval # Can happen if it's the first value!
             if idxs === nothing
@@ -118,11 +118,11 @@ end
                 vals[j] = u[k][idxs]
             end
         else
-            typeof(id) <: SensitivityInterpolation && error(SENSITIVITY_INTERP_MESSAGE)
+            id isa SensitivityInterpolation && error(SENSITIVITY_INTERP_MESSAGE)
             dt = t[i] - t[i - 1]
             Θ = (tval - t[i - 1]) / dt
             idxs_internal = idxs
-            if typeof(id) <: HermiteInterpolation
+            if id isa HermiteInterpolation
                 vals[j] = interpolant(Θ, id, dt, u[i - 1], u[i], du[i - 1], du[i],
                     idxs_internal, deriv)
             else
@@ -143,7 +143,7 @@ times t (sorted), with values u and derivatives ks
     continuity::Symbol = :left) where {I, D}
     t = id.t
     u = id.u
-    typeof(id) <: HermiteInterpolation && (du = id.du)
+    id isa HermiteInterpolation && (du = id.du)
     tdir = sign(t[end] - t[1])
     idx = sortperm(tvals, rev = tdir < 0)
     i = 2 # Start the search thinking it's between t[1] and t[2]
@@ -156,7 +156,7 @@ times t (sorted), with values u and derivatives ks
     for j in idx
         tval = tvals[j]
         i = searchsortedfirst(@view(t[i:end]), tval, rev = tdir < 0) + i - 1 # It's in the interval t[i-1] to t[i]
-        avoid_constant_ends = deriv != Val{0} #|| typeof(tval) <: ForwardDiff.Dual
+        avoid_constant_ends = deriv != Val{0} #|| tval isa ForwardDiff.Dual
         avoid_constant_ends && i == 1 && (i += 1)
         if !avoid_constant_ends && t[i - 1] == tval # Can happen if it's the first value!
             if idxs === nothing
@@ -173,19 +173,19 @@ times t (sorted), with values u and derivatives ks
                 vals[j] = u[k][idxs]
             end
         else
-            typeof(id) <: SensitivityInterpolation && error(SENSITIVITY_INTERP_MESSAGE)
+            id isa SensitivityInterpolation && error(SENSITIVITY_INTERP_MESSAGE)
             dt = t[i] - t[i - 1]
             Θ = (tval - t[i - 1]) / dt
             idxs_internal = idxs
             if eltype(u) <: Union{AbstractArray, ArrayPartition}
-                if typeof(id) <: HermiteInterpolation
+                if id isa HermiteInterpolation
                     interpolant!(vals[j], Θ, id, dt, u[i - 1], u[i], du[i - 1], du[i],
                         idxs_internal, deriv)
                 else
                     interpolant!(vals[j], Θ, id, dt, u[i - 1], u[i], idxs_internal, deriv)
                 end
             else
-                if typeof(id) <: HermiteInterpolation
+                if id isa HermiteInterpolation
                     vals[j] = interpolant(Θ, id, dt, u[i - 1], u[i], du[i - 1], du[i],
                         idxs_internal, deriv)
                 else
@@ -206,7 +206,7 @@ times t (sorted), with values u and derivatives ks
     continuity::Symbol = :left) where {I, D}
     t = id.t
     u = id.u
-    typeof(id) <: HermiteInterpolation && (du = id.du)
+    id isa HermiteInterpolation && (du = id.du)
     tdir = sign(t[end] - t[1])
     t[end] == t[1] && tval != t[end] &&
         error("Solution interpolation cannot extrapolate from a single timepoint. Either solve on a longer timespan or use the local extrapolation from the integrator interface.")
@@ -215,7 +215,7 @@ times t (sorted), with values u and derivatives ks
     tdir * tval < tdir * t[1] &&
         error("Solution interpolation cannot extrapolate before the first timepoint. Either start solving earlier or use the local extrapolation from the integrator interface.")
     @inbounds i = searchsortedfirst(t, tval, rev = tdir < 0) # It's in the interval t[i-1] to t[i]
-    avoid_constant_ends = deriv != Val{0} #|| typeof(tval) <: ForwardDiff.Dual
+    avoid_constant_ends = deriv != Val{0} #|| tval isa ForwardDiff.Dual
     avoid_constant_ends && i == 1 && (i += 1)
     if !avoid_constant_ends && t[i] == tval
         lasti = lastindex(t)
@@ -232,11 +232,11 @@ times t (sorted), with values u and derivatives ks
             val = u[i - 1][idxs]
         end
     else
-        typeof(id) <: SensitivityInterpolation && error(SENSITIVITY_INTERP_MESSAGE)
+        id isa SensitivityInterpolation && error(SENSITIVITY_INTERP_MESSAGE)
         dt = t[i] - t[i - 1]
         Θ = (tval - t[i - 1]) / dt
         idxs_internal = idxs
-        if typeof(id) <: HermiteInterpolation
+        if id isa HermiteInterpolation
             val = interpolant(Θ, id, dt, u[i - 1], u[i], du[i - 1], du[i], idxs_internal,
                 deriv)
         else
@@ -256,7 +256,7 @@ times t (sorted), with values u and derivatives ks
     continuity::Symbol = :left) where {I, D}
     t = id.t
     u = id.u
-    typeof(id) <: HermiteInterpolation && (du = id.du)
+    id isa HermiteInterpolation && (du = id.du)
     tdir = sign(t[end] - t[1])
     t[end] == t[1] && tval != t[end] &&
         error("Solution interpolation cannot extrapolate from a single timepoint. Either solve on a longer timespan or use the local extrapolation from the integrator interface.")
@@ -265,7 +265,7 @@ times t (sorted), with values u and derivatives ks
     tdir * tval < tdir * t[1] &&
         error("Solution interpolation cannot extrapolate before the first timepoint. Either start solving earlier or use the local extrapolation from the integrator interface.")
     @inbounds i = searchsortedfirst(t, tval, rev = tdir < 0) # It's in the interval t[i-1] to t[i]
-    avoid_constant_ends = deriv != Val{0} #|| typeof(tval) <: ForwardDiff.Dual
+    avoid_constant_ends = deriv != Val{0} #|| tval isa ForwardDiff.Dual
     avoid_constant_ends && i == 1 && (i += 1)
     if !avoid_constant_ends && t[i] == tval
         lasti = lastindex(t)
@@ -282,11 +282,11 @@ times t (sorted), with values u and derivatives ks
             copy!(out, u[i - 1][idxs])
         end
     else
-        typeof(id) <: SensitivityInterpolation && error(SENSITIVITY_INTERP_MESSAGE)
+        id isa SensitivityInterpolation && error(SENSITIVITY_INTERP_MESSAGE)
         dt = t[i] - t[i - 1]
         Θ = (tval - t[i - 1]) / dt
         idxs_internal = idxs
-        if typeof(id) <: HermiteInterpolation
+        if id isa HermiteInterpolation
             interpolant!(out, Θ, id, dt, u[i - 1], u[i], du[i - 1], du[i], idxs_internal,
                 deriv)
         else
diff --git a/src/operators/diffeq_operator.jl b/src/operators/diffeq_operator.jl
index e508c7022..a97063b83 100644
--- a/src/operators/diffeq_operator.jl
+++ b/src/operators/diffeq_operator.jl
@@ -37,7 +37,7 @@ function (L::AffineDiffEqOperator)(u, p, t::Number)
     update_coefficients!(L, u, p, t)
     du = sum(A * u for A in L.As)
     for B in L.Bs
-        if typeof(B) <: Union{Number, AbstractArray}
+        if B isa Union{Number, AbstractArray}
             du .+= B
         else
             du .+= B(t)
@@ -58,7 +58,7 @@ function (L::AffineDiffEqOperator)(du, u, p, t::Number)
         du .+= du_cache
     end
     for B in L.Bs
-        if typeof(B) <: Union{Number, AbstractArray}
+        if B isa Union{Number, AbstractArray}
             du .+= B
         else
             B(du_cache, t)
diff --git a/src/problems/discrete_problems.jl b/src/problems/discrete_problems.jl
index b2f2a362c..4b06bdfa5 100644
--- a/src/problems/discrete_problems.jl
+++ b/src/problems/discrete_problems.jl
@@ -147,7 +147,7 @@ Define a discrete problem with the identity map.
 """
 function DiscreteProblem(u0::Union{AbstractArray, Number}, tspan::Tuple,
     p = NullParameters(); kwargs...)
-    iip = typeof(u0) <: AbstractArray
+    iip = u0 isa AbstractArray
     if iip
         f = DISCRETE_INPLACE_DEFAULT
     else
diff --git a/src/problems/problem_utils.jl b/src/problems/problem_utils.jl
index 936eab5b7..c72a451f3 100644
--- a/src/problems/problem_utils.jl
+++ b/src/problems/problem_utils.jl
@@ -23,7 +23,7 @@ function Base.summary(io::IO, prob::AbstractDEProblem)
         type_color, typeof(prob.u0),
         no_color, " and tType ",
         type_color,
-        typeof(prob.tspan) <: Function ?
+        prob.tspan isa Function ?
         "Unknown" : (prob.tspan === nothing ?
          "Nothing" : typeof(prob.tspan[1])),
         no_color, ". In-place: ",
diff --git a/src/remake.jl b/src/remake.jl
index 083f0fda6..a03c11ec5 100644
--- a/src/remake.jl
+++ b/src/remake.jl
@@ -403,7 +403,7 @@ end
 
 # overloaded in MTK to intercept symbolic remake
 function process_p_u0_symbolic(prob, p, u0)
-    if typeof(prob) <: Union{AbstractDEProblem, OptimizationProblem, NonlinearProblem}
+    if prob isa Union{AbstractDEProblem, OptimizationProblem, NonlinearProblem}
         throw(ArgumentError("Please load `ModelingToolkit.jl` in order to support symbolic remake."))
     else
         throw(ArgumentError("Symbolic remake for $(typeof(prob)) is currently not supported, consider opening an issue."))
diff --git a/src/scimlfunctions.jl b/src/scimlfunctions.jl
index 6afcbf108..55eae6b4b 100644
--- a/src/scimlfunctions.jl
+++ b/src/scimlfunctions.jl
@@ -2488,7 +2488,7 @@ function ODEFunction{iip, specialize}(f;
     sys = __has_sys(f) ? f.sys : nothing) where {iip,
     specialize,
 }
-    if mass_matrix === I && typeof(f) <: Tuple
+    if mass_matrix === I && f isa Tuple
         mass_matrix = ((I for i in 1:length(f))...,)
     end
 
@@ -2705,7 +2705,7 @@ end
     f1 = ODEFunction(f1)
     f2 = ODEFunction(f2)
 
-    if !(typeof(f1) <: AbstractSciMLOperator || typeof(f1.f) <: AbstractSciMLOperator) &&
+    if !(f1 isa AbstractSciMLOperator || f1.f isa AbstractSciMLOperator) &&
        isinplace(f1) != isinplace(f2)
         throw(NonconformingFunctionsError(["f2"]))
     end
@@ -2787,7 +2787,7 @@ SplitFunction(f::SplitFunction; kwargs...) = f
     jvp, vjp, jac_prototype, sparsity, Wfact,
     Wfact_t, paramjac, syms, indepsym, paramsyms,
     observed, colorvec, sys) where {iip}
-    f1 = typeof(f1) <: AbstractSciMLOperator ? f1 : ODEFunction(f1)
+    f1 = f1 isa AbstractSciMLOperator ? f1 : ODEFunction(f1)
     f2 = ODEFunction(f2)
 
     if isinplace(f1) != isinplace(f2)
@@ -3151,7 +3151,7 @@ SDEFunction(f::SDEFunction; kwargs...) = f
     jvp, vjp,
     jac_prototype, Wfact, Wfact_t, paramjac, observed,
     syms, indepsym, paramsyms, colorvec, sys)
-    f1 = typeof(f1) <: AbstractSciMLOperator ? f1 : SDEFunction(f1)
+    f1 = f1 isa AbstractSciMLOperator ? f1 : SDEFunction(f1)
     f2 = SDEFunction(f2)
     SplitFunction{isinplace(f2), typeof(f1), typeof(f2), typeof(g), typeof(mass_matrix),
         typeof(cache), typeof(analytic), typeof(tgrad), typeof(jac), typeof(jvp),
@@ -3237,7 +3237,7 @@ SplitSDEFunction(f::SplitSDEFunction; kwargs...) = f
     jac_prototype, Wfact, Wfact_t, paramjac,
     syms, indepsym, paramsyms, observed, colorvec,
     sys)
-    f1 = typeof(f1) <: AbstractSciMLOperator ? f1 : SDEFunction(f1)
+    f1 = f1 isa AbstractSciMLOperator ? f1 : SDEFunction(f1)
     f2 = SDEFunction(f2)
     DynamicalSDEFunction{isinplace(f2), FullSpecialize, typeof(f1), typeof(f2), typeof(g),
         typeof(mass_matrix),
@@ -3603,7 +3603,7 @@ DDEFunction(f::DDEFunction; kwargs...) = f
     paramjac,
     syms, indepsym, paramsyms, observed,
     colorvec) where {iip}
-    f1 = typeof(f1) <: AbstractSciMLOperator ? f1 : DDEFunction(f1)
+    f1 = f1 isa AbstractSciMLOperator ? f1 : DDEFunction(f1)
     f2 = DDEFunction(f2)
     DynamicalDDEFunction{isinplace(f2), FullSpecialize, typeof(f1), typeof(f2),
         typeof(mass_matrix),
@@ -3808,7 +3808,7 @@ function NonlinearFunction{iip, specialize}(f;
     resid_prototype = __has_resid_prototype(f) ? f.resid_prototype : nothing) where {
     iip, specialize}
 
-    if mass_matrix === I && typeof(f) <: Tuple
+    if mass_matrix === I && f isa Tuple
         mass_matrix = ((I for i in 1:length(f))...,)
     end
 
@@ -3995,7 +3995,7 @@ function BVPFunction{iip, specialize, twopoint}(f, bc;
     colorvec = __has_colorvec(f) ? f.colorvec : nothing,
     bccolorvec = __has_colorvec(bc) ? bc.colorvec : nothing,
     sys = __has_sys(f) ? f.sys : nothing) where {iip, specialize, twopoint}
-    if mass_matrix === I && typeof(f) <: Tuple
+    if mass_matrix === I && f isa Tuple
         mass_matrix = ((I for i in 1:length(f))...,)
     end
 
diff --git a/src/solutions/ode_solutions.jl b/src/solutions/ode_solutions.jl
index 77f13cc39..b3ba74fe6 100644
--- a/src/solutions/ode_solutions.jl
+++ b/src/solutions/ode_solutions.jl
@@ -261,7 +261,7 @@ function build_solution(prob::Union{AbstractODEProblem, AbstractDDEProblem},
         N = length((size(prob.u0)..., length(u)))
     end
 
-    if typeof(prob.f) <: Tuple
+    if prob.f isa Tuple
         f = prob.f[1]
     else
         f = prob.f
@@ -419,9 +419,9 @@ end
 function sensitivity_solution(sol::ODESolution, u, t)
     T = eltype(eltype(u))
     N = length((size(sol.prob.u0)..., length(u)))
-    interp = if typeof(sol.interp) <: LinearInterpolation
+    interp = if sol.interp isa LinearInterpolation
         LinearInterpolation(t, u)
-    elseif typeof(sol.interp) <: ConstantInterpolation
+    elseif sol.interp isa ConstantInterpolation
         ConstantInterpolation(t, u)
     else
         SensitivityInterpolation(t, u)
diff --git a/src/solutions/rode_solutions.jl b/src/solutions/rode_solutions.jl
index 6597ac185..350c3c823 100644
--- a/src/solutions/rode_solutions.jl
+++ b/src/solutions/rode_solutions.jl
@@ -80,7 +80,7 @@ function build_solution(prob::Union{AbstractRODEProblem, AbstractSDDEProblem},
     T = eltype(eltype(u))
     N = length((size(prob.u0)..., length(u)))
 
-    if typeof(prob.f) <: Tuple
+    if prob.f isa Tuple
         f = prob.f[1]
     else
         f = prob.f
@@ -134,7 +134,7 @@ end
 
 function calculate_solution_errors!(sol::AbstractRODESolution; fill_uanalytic = true,
     timeseries_errors = true, dense_errors = true)
-    if typeof(sol.prob.f) <: Tuple
+    if sol.prob.f isa Tuple
         f = sol.prob.f[1]
     else
         f = sol.prob.f
@@ -232,9 +232,9 @@ end
 function sensitivity_solution(sol::AbstractRODESolution, u, t)
     T = eltype(eltype(u))
     N = length((size(sol.prob.u0)..., length(u)))
-    interp = if typeof(sol.interp) <: LinearInterpolation
+    interp = if sol.interp isa LinearInterpolation
         LinearInterpolation(t, u)
-    elseif typeof(sol.interp) <: ConstantInterpolation
+    elseif sol.interp isa ConstantInterpolation
         ConstantInterpolation(t, u)
     else
         SensitivityInterpolation(t, u)
diff --git a/src/solutions/solution_interface.jl b/src/solutions/solution_interface.jl
index 900fecfea..ed0f5f88c 100644
--- a/src/solutions/solution_interface.jl
+++ b/src/solutions/solution_interface.jl
@@ -233,13 +233,13 @@ DEFAULT_PLOT_FUNC(x, y, z) = (x, y, z) # For v0.5.2 bug
 @recipe function f(sol::AbstractTimeseriesSolution;
     plot_analytic = false,
     denseplot = (sol.dense ||
-                 typeof(sol.prob) <: AbstractDiscreteProblem) &&
-                !(typeof(sol) <: AbstractRODESolution) &&
+                 sol.prob isa AbstractDiscreteProblem) &&
+                !(sol isa AbstractRODESolution) &&
                 !(hasfield(typeof(sol), :interp) &&
-                  typeof(sol.interp) <: SensitivityInterpolation),
+                  sol.interp isa SensitivityInterpolation),
     plotdensity = min(Int(1e5),
         sol.tslocation == 0 ?
-        (typeof(sol.prob) <: AbstractDiscreteProblem ?
+        (sol.prob isa AbstractDiscreteProblem ?
          max(1000, 100 * length(sol)) :
          max(1000, 10 * length(sol))) :
         1000 * sol.tslocation),
@@ -271,7 +271,7 @@ DEFAULT_PLOT_FUNC(x, y, z) = (x, y, z) # For v0.5.2 bug
     seriestype --> :path
 
     # Special case labels when idxs = (:x,:y,:z) or (:x) or [:x,:y] ...
-    if typeof(idxs) <: Tuple && (issymbollike(idxs[1]) && issymbollike(idxs[2]))
+    if idxs isa Tuple && (issymbollike(idxs[1]) && issymbollike(idxs[2]))
         val = issymbollike(int_vars[1][2]) ? String(Symbol(int_vars[1][2])) :
               strs[int_vars[1][2]]
         xguide --> val
@@ -330,7 +330,7 @@ DEFAULT_PLOT_FUNC(x, y, z) = (x, y, z) # For v0.5.2 bug
     # Analytical solutions do not save enough information to have a good idea
     # of the axis ahead of time
     # Only set axis for animations
-    if sol.tslocation != 0 && !(typeof(sol) <: AbstractAnalyticalSolution)
+    if sol.tslocation != 0 && !(sol isa AbstractAnalyticalSolution)
         if all(getindex.(int_vars, 1) .== DEFAULT_PLOT_FUNC)
             mins = minimum(sol[int_vars[1][3], :])
             maxs = maximum(sol[int_vars[1][3], :])
@@ -381,9 +381,9 @@ function diffeq_to_arrays(sol, plot_analytic, denseplot, plotdensity, tspan, axi
 
     if denseplot
         # Generate the points from the plot from dense function
-        if tspan === nothing && !(typeof(sol) <: AbstractAnalyticalSolution)
+        if tspan === nothing && !(sol isa AbstractAnalyticalSolution)
             plott = collect(densetspacer(sol.t[start_idx], sol.t[end_idx], plotdensity))
-        elseif typeof(sol) <: AbstractAnalyticalSolution
+        elseif sol isa AbstractAnalyticalSolution
             tspan = sol.prob.tspan
             plott = collect(densetspacer(tspan[1], tspan[end], plotdensity))
         else
@@ -391,7 +391,7 @@ function diffeq_to_arrays(sol, plot_analytic, denseplot, plotdensity, tspan, axi
         end
         plot_timeseries = sol(plott)
         if plot_analytic
-            if typeof(sol.prob.f) <: Tuple
+            if sol.prob.f isa Tuple
                 plot_analytic_timeseries = [sol.prob.f[1].analytic(sol.prob.u0, sol.prob.p,
                     t) for t in plott]
             else
@@ -442,7 +442,7 @@ function interpret_vars(vars, sol, syms)
         # Do syms conversion
         tmp_vars = []
         for var in vars
-            if typeof(var) <: Union{Tuple, AbstractArray} #eltype(var) <: Symbol # Some kind of iterable
+            if var isa Union{Tuple, AbstractArray} #eltype(var) <: Symbol # Some kind of iterable
                 tmp = []
                 for x in var
                     if issymbollike(x)
@@ -454,7 +454,7 @@ function interpret_vars(vars, sol, syms)
                         push!(tmp, x)
                     end
                 end
-                if typeof(var) <: Tuple
+                if var isa Tuple
                     var_int = tuple(tmp...)
                 else
                     var_int = tmp
@@ -477,7 +477,7 @@ function interpret_vars(vars, sol, syms)
             end
             push!(tmp_vars, var_int)
         end
-        if typeof(vars) <: Tuple
+        if vars isa Tuple
             vars = tuple(tmp_vars...)
         else
             vars = tmp_vars
@@ -486,23 +486,23 @@ function interpret_vars(vars, sol, syms)
 
     if vars === nothing
         # Default: plot all timeseries
-        if typeof(sol[1]) <: Union{Tuple, AbstractArray}
+        if sol[1] isa Union{Tuple, AbstractArray}
             vars = collect((DEFAULT_PLOT_FUNC, 0, i) for i in plot_indices(sol[1]))
         else
             vars = [(DEFAULT_PLOT_FUNC, 0, 1)]
         end
     end
 
-    if typeof(vars) <: Base.Integer
+    if vars isa Base.Integer
         vars = [(DEFAULT_PLOT_FUNC, 0, vars)]
     end
 
-    if typeof(vars) <: AbstractArray
+    if vars isa AbstractArray
         # If list given, its elements should be tuples, or we assume x = time
         tmp = Tuple[]
         for x in vars
-            if typeof(x) <: Tuple
-                if typeof(x[1]) <: Int
+            if x isa Tuple
+                if x[1] isa Int
                     push!(tmp, tuple(DEFAULT_PLOT_FUNC, x...))
                 else
                     push!(tmp, x)
@@ -514,10 +514,10 @@ function interpret_vars(vars, sol, syms)
         vars = tmp
     end
 
-    if typeof(vars) <: Tuple
+    if vars isa Tuple
         # If tuple given...
-        if typeof(vars[end - 1]) <: AbstractArray
-            if typeof(vars[end]) <: AbstractArray
+        if vars[end - 1] isa AbstractArray
+            if vars[end] isa AbstractArray
                 # If both axes are lists we zip (will fail if different lengths)
                 vars = collect(zip([DEFAULT_PLOT_FUNC for i in eachindex(vars[end - 1])],
                     vars[end - 1], vars[end]))
@@ -526,12 +526,12 @@ function interpret_vars(vars, sol, syms)
                 vars = [(DEFAULT_PLOT_FUNC, x, vars[end]) for x in vars[end - 1]]
             end
         else
-            if typeof(vars[2]) <: AbstractArray
+            if vars[2] isa AbstractArray
                 # Just the y axis is a list
                 vars = [(DEFAULT_PLOT_FUNC, vars[end - 1], y) for y in vars[end]]
             else
                 # Both axes are numbers
-                if typeof(vars[1]) <: Int || issymbollike(vars[1])
+                if vars[1] isa Int || issymbollike(vars[1])
                     vars = [tuple(DEFAULT_PLOT_FUNC, vars...)]
                 else
                     vars = [vars]
@@ -564,7 +564,7 @@ function add_labels!(labels, x, dims, sol, strs)
     lys[end] = chop(lys[end]) # Take off the last comma
     if !issymbollike(x[2]) && x[2] == 0 && dims == 3
         # if there are no dependence in syms, then we add "(t)"
-        if strs !== nothing && (typeof(x[3]) <: Int && endswith(strs[x[3]], r"(.*)")) ||
+        if strs !== nothing && (x[3] isa Int && endswith(strs[x[3]], r"(.*)")) ||
            (issymbollike(x[3]) && endswith(string(x[3]), r"(.*)"))
             tmp_lab = "$(lys...)"
         else
@@ -627,7 +627,7 @@ function u_n(timeseries::AbstractArray, n::Int, sol, plott, plot_timeseries)
     # Returns the nth variable from the timeseries, t if n == 0
     if n == 0
         return plott
-    elseif n == 1 && !(typeof(sol[1]) <: Union{AbstractArray, ArrayPartition})
+    elseif n == 1 && !(sol[1] isa Union{AbstractArray, ArrayPartition})
         return timeseries
     else
         tmp = Vector{eltype(sol[1])}(undef, length(plot_timeseries))