From 8ee5ca4db3a11430f73e098f33fe03919436c664 Mon Sep 17 00:00:00 2001
From: Avik Pal <avikpal@mit.edu>
Date: Thu, 10 Aug 2023 16:19:14 -0400
Subject: [PATCH 01/16] TwoPointBVPFunction restructure

---
 Project.toml                 |  2 +-
 src/problems/bvp_problems.jl | 10 +++++++---
 2 files changed, 8 insertions(+), 4 deletions(-)

diff --git a/Project.toml b/Project.toml
index a5efdc86d..ebf74000c 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,7 +1,7 @@
 name = "SciMLBase"
 uuid = "0bca4576-84f4-4d90-8ffe-ffa030f20462"
 authors = ["Chris Rackauckas <accounts@chrisrackauckas.com> and contributors"]
-version = "1.94.0"
+version = "2.0.0"
 
 [deps]
 ADTypes = "47edcb42-4c32-4615-8424-f2b9edc5f35b"
diff --git a/src/problems/bvp_problems.jl b/src/problems/bvp_problems.jl
index aa5d887b4..84c9715a4 100644
--- a/src/problems/bvp_problems.jl
+++ b/src/problems/bvp_problems.jl
@@ -133,9 +133,13 @@ $(TYPEDEF)
 struct TwoPointBVPFunction{bF}
     bc::bF
 end
-TwoPointBVPFunction(; bc = error("No argument bc")) = TwoPointBVPFunction(bc)
-(f::TwoPointBVPFunction)(residual, ua, ub, p) = f.bc(residual, ua, ub, p)
-(f::TwoPointBVPFunction)(residual, u, p) = f.bc(residual, u[1], u[end], p)
+TwoPointBVPFunction(; bc = error("No argument `bc`")) = TwoPointBVPFunction(bc)
+function (f::TwoPointBVPFunction)(residuala, residualb, ua, ub, p)
+    return f.bc(residuala, residualb, ua, ub, p)
+end
+function (f::TwoPointBVPFunction)(residual::Tuple, u, p)
+    return f(residual[1], residual[2], u[1], u[end], p)
+end
 
 """
 $(TYPEDEF)

From 75f258e34e2bc3a484ed772b37b18460f88e9b3a Mon Sep 17 00:00:00 2001
From: Avik Pal <avikpal@mit.edu>
Date: Tue, 5 Sep 2023 16:53:31 -0400
Subject: [PATCH 02/16] This is breaking

---
 Project.toml | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/Project.toml b/Project.toml
index a14d73463..951f6b7ca 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,7 +1,7 @@
 name = "SciMLBase"
 uuid = "0bca4576-84f4-4d90-8ffe-ffa030f20462"
 authors = ["Chris Rackauckas <accounts@chrisrackauckas.com> and contributors"]
-version = "1.96.1"
+version = "2.0.0"
 
 [deps]
 ADTypes = "47edcb42-4c32-4615-8424-f2b9edc5f35b"

From 7bfe639c3ab701b3a816fa9cbc958d2aad6502f0 Mon Sep 17 00:00:00 2001
From: Avik Pal <avikpal@mit.edu>
Date: Tue, 5 Sep 2023 17:02:49 -0400
Subject: [PATCH 03/16] Extend nonconforming to `bc`

---
 src/scimlfunctions.jl | 36 ++++++++++--------------------------
 1 file changed, 10 insertions(+), 26 deletions(-)

diff --git a/src/scimlfunctions.jl b/src/scimlfunctions.jl
index 102e0c041..d31677fef 100644
--- a/src/scimlfunctions.jl
+++ b/src/scimlfunctions.jl
@@ -2233,11 +2233,7 @@ For more details on this argument, see the ODEFunction documentation.
 The fields of the BVPFunction type directly match the names of the inputs.
 """
 struct BVPFunction{iip, specialize, F, BF, TMM, Ta, Tt, TJ, BCTJ, JVP, VJP, JP,
-    BCJP, SP, TW, TWt,
-    TPJ,
-    S, S2, S3, O, TCV, BCTCV,
-    SYS} <:
-       AbstractBVPFunction{iip}
+    BCJP, SP, TW, TWt, TPJ, S, S2, S3, O, TCV, BCTCV, SYS} <: AbstractBVPFunction{iip}
     f::F
     bc::BF
     mass_matrix::TMM
@@ -3815,31 +3811,23 @@ function OptimizationFunction{iip}(f, adtype::AbstractADType = NoAD();
 end
 
 function BVPFunction{iip, specialize}(f, bc;
-    mass_matrix = __has_mass_matrix(f) ? f.mass_matrix :
-                  I,
+    mass_matrix = __has_mass_matrix(f) ? f.mass_matrix : I,
     analytic = __has_analytic(f) ? f.analytic : nothing,
     tgrad = __has_tgrad(f) ? f.tgrad : nothing,
     jac = __has_jac(f) ? f.jac : nothing,
     bcjac = __has_jac(bc) ? bc.jac : nothing,
     jvp = __has_jvp(f) ? f.jvp : nothing,
     vjp = __has_vjp(f) ? f.vjp : nothing,
-    jac_prototype = __has_jac_prototype(f) ?
-                    f.jac_prototype :
-                    nothing,
-    bcjac_prototype = __has_jac_prototype(bc) ?
-                      bc.jac_prototype :
-                      nothing,
-    sparsity = __has_sparsity(f) ? f.sparsity :
-               jac_prototype,
+    jac_prototype = __has_jac_prototype(f) ? f.jac_prototype : nothing,
+    bcjac_prototype = __has_jac_prototype(bc) ? bc.jac_prototype : nothing,
+    sparsity = __has_sparsity(f) ? f.sparsity :            jac_prototype,
     Wfact = __has_Wfact(f) ? f.Wfact : nothing,
     Wfact_t = __has_Wfact_t(f) ? f.Wfact_t : nothing,
     paramjac = __has_paramjac(f) ? f.paramjac : nothing,
     syms = __has_syms(f) ? f.syms : nothing,
     indepsym = __has_indepsym(f) ? f.indepsym : nothing,
-    paramsyms = __has_paramsyms(f) ? f.paramsyms :
-                nothing,
-    observed = __has_observed(f) ? f.observed :
-               DEFAULT_OBSERVED,
+    paramsyms = __has_paramsyms(f) ? f.paramsyms : nothing,
+    observed = __has_observed(f) ? f.observed : DEFAULT_OBSERVED,
     colorvec = __has_colorvec(f) ? f.colorvec : nothing,
     bccolorvec = __has_colorvec(bc) ? bc.colorvec : nothing,
     sys = __has_sys(f) ? f.sys : nothing) where {iip, specialize}
@@ -3892,17 +3880,13 @@ function BVPFunction{iip, specialize}(f, bc;
     Wfact_tiip = Wfact_t !== nothing ? isinplace(Wfact_t, 5, "Wfact_t", iip) : iip
     paramjaciip = paramjac !== nothing ? isinplace(paramjac, 4, "paramjac", iip) : iip
 
-    nonconforming = (jaciip,
-        tgradiip,
-        jvpiip,
-        vjpiip,
-        Wfactiip,
-        Wfact_tiip,
+    nonconforming = (bciip, jaciip, tgradiip, jvpiip, vjpiip, Wfactiip, Wfact_tiip,
         paramjaciip) .!= iip
     bc_nonconforming = bcjaciip .!= bciip
     if any(nonconforming)
         nonconforming = findall(nonconforming)
-        functions = ["jac", "bcjac", "tgrad", "jvp", "vjp", "Wfact", "Wfact_t", "paramjac"][nonconforming]
+        functions = ["bc", "jac", "bcjac", "tgrad", "jvp", "vjp", "Wfact", "Wfact_t",
+            "paramjac"][nonconforming]
         throw(NonconformingFunctionsError(functions))
     end
 

From e47627367861475f3625d40b91a3f38b106bb2d3 Mon Sep 17 00:00:00 2001
From: Avik Pal <avikpal@mit.edu>
Date: Tue, 5 Sep 2023 18:11:55 -0400
Subject: [PATCH 04/16] Don't dump BVPFunction Information

---
 src/problems/bvp_problems.jl | 42 +++++++++++++++---------------------
 src/scimlfunctions.jl        | 18 ++++++++++------
 2 files changed, 28 insertions(+), 32 deletions(-)

diff --git a/src/problems/bvp_problems.jl b/src/problems/bvp_problems.jl
index a25ea8f74..d9787f9bc 100644
--- a/src/problems/bvp_problems.jl
+++ b/src/problems/bvp_problems.jl
@@ -89,15 +89,12 @@ struct BVProblem{uType, tType, isinplace, P, F, BF, PT, K} <:
     kwargs::K
 
     @add_kwonly function BVProblem{iip}(f::AbstractBVPFunction{iip}, bc, u0, tspan,
-        p = NullParameters(),
-        problem_type = StandardBVProblem();
-        kwargs...) where {iip}
+        p = NullParameters(), problem_type = StandardBVProblem(); kwargs...) where {iip}
         _tspan = promote_tspan(tspan)
         warn_paramtype(p)
-        new{typeof(u0), typeof(_tspan), iip, typeof(p),
-            typeof(f.f), typeof(bc),
-            typeof(problem_type), typeof(kwargs)}(f.f, bc, u0, _tspan, p,
-            problem_type, kwargs)
+        return new{typeof(u0), typeof(_tspan), iip, typeof(p), typeof(f), typeof(bc),
+            typeof(problem_type), typeof(kwargs)}(f, bc, u0, _tspan, p, problem_type,
+            kwargs)
     end
 
     function BVProblem{iip}(f, bc, u0, tspan, p = NullParameters(); kwargs...) where {iip}
@@ -108,11 +105,12 @@ end
 TruncatedStacktraces.@truncate_stacktrace BVProblem 3 1 2
 
 function BVProblem(f, bc, u0, tspan, p = NullParameters(); kwargs...)
-    BVProblem(BVPFunction(f, bc), u0, tspan, p; kwargs...)
+    iip = isinplace(f, 4)
+    return BVProblem{iip}(BVPFunction{iip}(f, bc), u0, tspan, p; kwargs...)
 end
 
 function BVProblem(f::AbstractBVPFunction, u0, tspan, p = NullParameters(); kwargs...)
-    BVProblem{isinplace(f)}(f.f, f.bc, u0, tspan, p; kwargs...)
+    return BVProblem{isinplace(f)}(f, f.bc, u0, tspan, p; kwargs...)
 end
 
 """
@@ -132,31 +130,25 @@ end
 """
 $(TYPEDEF)
 """
-struct TwoPointBVProblem{iip} end
+struct TwoPointBVProblem end
+
 function TwoPointBVProblem(f, bc, u0, tspan, p = NullParameters(); kwargs...)
-    iip = isinplace(f, 4)
-    TwoPointBVProblem{iip}(f, bc, u0, tspan, p; kwargs...)
+    return TwoPointBVProblem(BVPFunction(f, bc), u0, tspan, p; kwargs...)
 end
-function TwoPointBVProblem{iip}(f, bc, u0, tspan, p = NullParameters();
-    kwargs...) where {iip}
-    BVProblem{iip}(f, TwoPointBVPFunction(bc), u0, tspan, p; kwargs...)
+function TwoPointBVProblem(f::AbstractBVPFunction, u0, tspan, p = NullParameters();
+    kwargs...)
+    iip = isinplace(f)
+    return BVProblem{iip}(f, f.bc, u0, tspan, p, TwoPointBVProblem(); kwargs...)
 end
 
 # Allow previous timeseries solution
-function TwoPointBVProblem(f::AbstractODEFunction,
-    bc,
-    sol::T,
-    tspan::Tuple,
+function TwoPointBVProblem(f::AbstractODEFunction, bc, sol::T, tspan::Tuple,
     p = NullParameters()) where {T <: AbstractTimeseriesSolution}
     TwoPointBVProblem(f, bc, sol.u, tspan, p)
 end
 # Allow initial guess function for the initial guess
-function TwoPointBVProblem(f::AbstractODEFunction,
-    bc,
-    initialGuess,
-    tspan::AbstractVector,
-    p = NullParameters();
-    kwargs...)
+function TwoPointBVProblem(f::AbstractODEFunction, bc, initialGuess, tspan::AbstractVector,
+    p = NullParameters(); kwargs...)
     u0 = [initialGuess(i) for i in tspan]
     TwoPointBVProblem(f, bc, u0, (tspan[1], tspan[end]), p)
 end
diff --git a/src/scimlfunctions.jl b/src/scimlfunctions.jl
index d31677fef..6cd7a6b4b 100644
--- a/src/scimlfunctions.jl
+++ b/src/scimlfunctions.jl
@@ -2233,7 +2233,7 @@ For more details on this argument, see the ODEFunction documentation.
 The fields of the BVPFunction type directly match the names of the inputs.
 """
 struct BVPFunction{iip, specialize, F, BF, TMM, Ta, Tt, TJ, BCTJ, JVP, VJP, JP,
-    BCJP, SP, TW, TWt, TPJ, S, S2, S3, O, TCV, BCTCV, SYS} <: AbstractBVPFunction{iip}
+    BCJP, BCRP, SP, TW, TWt, TPJ, S, S2, S3, O, TCV, BCTCV, SYS} <: AbstractBVPFunction{iip}
     f::F
     bc::BF
     mass_matrix::TMM
@@ -2245,6 +2245,7 @@ struct BVPFunction{iip, specialize, F, BF, TMM, Ta, Tt, TJ, BCTJ, JVP, VJP, JP,
     vjp::VJP
     jac_prototype::JP
     bcjac_prototype::BCJP
+    bcresid_prototype::BCRP
     sparsity::SP
     Wfact::TW
     Wfact_t::TWt
@@ -3820,7 +3821,8 @@ function BVPFunction{iip, specialize}(f, bc;
     vjp = __has_vjp(f) ? f.vjp : nothing,
     jac_prototype = __has_jac_prototype(f) ? f.jac_prototype : nothing,
     bcjac_prototype = __has_jac_prototype(bc) ? bc.jac_prototype : nothing,
-    sparsity = __has_sparsity(f) ? f.sparsity :            jac_prototype,
+    bcresid_prototype = nothing,
+    sparsity = __has_sparsity(f) ? f.sparsity : jac_prototype,
     Wfact = __has_Wfact(f) ? f.Wfact : nothing,
     Wfact_t = __has_Wfact_t(f) ? f.Wfact_t : nothing,
     paramjac = __has_paramjac(f) ? f.paramjac : nothing,
@@ -3831,6 +3833,7 @@ function BVPFunction{iip, specialize}(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}
+    # TODO: Add checks for bcresid_prototype for TwoPointBVPFunction
     if mass_matrix === I && typeof(f) <: Tuple
         mass_matrix = ((I for i in 1:length(f))...,)
     end
@@ -3898,14 +3901,14 @@ function BVPFunction{iip, specialize}(f, bc;
 
     if specialize === NoSpecialize
         BVPFunction{iip, specialize, Any, Any, Any, Any, Any,
-            Any, Any, Any, Any, Any, Any, Any, Any, Any,
+            Any, Any, Any, Any, Any, Any, Any, Any, Any, Any,
             Any, typeof(syms), typeof(indepsym), typeof(paramsyms),
             Any, typeof(_colorvec), typeof(_bccolorvec), Any}(f, bc, mass_matrix,
             analytic,
             tgrad,
             jac, bcjac, jvp, vjp,
             jac_prototype,
-            bcjac_prototype,
+            bcjac_prototype, bcresid_prototype,
             sparsity, Wfact,
             Wfact_t,
             paramjac, syms,
@@ -3916,11 +3919,12 @@ function BVPFunction{iip, specialize}(f, bc;
         BVPFunction{iip, FunctionWrapperSpecialize,
             typeof(f), typeof(bc), typeof(mass_matrix), typeof(analytic), typeof(tgrad),
             typeof(jac), typeof(jvp), typeof(vjp), typeof(jac_prototype),
+            typeof(bcjac_prototype), typeof(bcresid_prototype),
             typeof(sparsity), typeof(Wfact), typeof(Wfact_t), typeof(paramjac),
             typeof(syms), typeof(indepsym), typeof(paramsyms), typeof(observed),
             typeof(_colorvec), typeof(_bccolorvec),
             typeof(sys)}(f, bc, mass_matrix, analytic, tgrad, jac, bcjac,
-            jvp, vjp, jac_prototype, bcjac_prototype, sparsity, Wfact,
+            jvp, vjp, jac_prototype, bcjac_prototype, bcresid_prototype, sparsity, Wfact,
             Wfact_t, paramjac, syms, indepsym, paramsyms,
             observed, _colorvec, _bccolorvec, sys)
     else
@@ -3928,13 +3932,13 @@ function BVPFunction{iip, specialize}(f, bc;
             typeof(analytic),
             typeof(tgrad),
             typeof(jac), typeof(bcjac), typeof(jvp), typeof(vjp), typeof(jac_prototype),
-            typeof(bcjac_prototype),
+            typeof(bcjac_prototype), typeof(bcresid_prototype),
             typeof(sparsity), typeof(Wfact), typeof(Wfact_t),
             typeof(paramjac), typeof(syms), typeof(indepsym), typeof(paramsyms),
             typeof(observed),
             typeof(_colorvec), typeof(_bccolorvec), typeof(sys)}(f, bc, mass_matrix, analytic,
             tgrad, jac, bcjac, jvp, vjp,
-            jac_prototype, bcjac_prototype, sparsity,
+            jac_prototype, bcjac_prototype, bcresid_prototype, sparsity,
             Wfact, Wfact_t, paramjac,
             syms, indepsym, paramsyms, observed,
             _colorvec, _bccolorvec, sys)

From 2b73aefd267af7f213d8fc2f5ef72e3083d213ae Mon Sep 17 00:00:00 2001
From: Avik Pal <avikpal@mit.edu>
Date: Wed, 6 Sep 2023 14:16:21 -0400
Subject: [PATCH 05/16] Fix TwoPointBVP Dispatches

---
 src/problems/bvp_problems.jl | 81 +++++++++++++++++++++---------------
 src/scimlfunctions.jl        | 71 ++++++++++++-------------------
 2 files changed, 73 insertions(+), 79 deletions(-)

diff --git a/src/problems/bvp_problems.jl b/src/problems/bvp_problems.jl
index d9787f9bc..b46b0479c 100644
--- a/src/problems/bvp_problems.jl
+++ b/src/problems/bvp_problems.jl
@@ -3,6 +3,11 @@ $(TYPEDEF)
 """
 struct StandardBVProblem end
 
+"""
+$(TYPEDEF)
+"""
+struct TwoPointBVProblem end
+
 @doc doc"""
 
 Defines an BVP problem.
@@ -17,7 +22,7 @@ condition ``u_0`` which define an ODE:
 \frac{du}{dt} = f(u,p,t)
 ```
 
-along with an implicit function `bc!` which defines the residual equation, where
+along with an implicit function `bc` which defines the residual equation, where
 
 ```math
 bc(u,p,t) = 0
@@ -36,22 +41,27 @@ u(t_f) = b
 ### Constructors
 
 ```julia
-TwoPointBVProblem{isinplace}(f,bc!,u0,tspan,p=NullParameters();kwargs...)
-BVProblem{isinplace}(f,bc!,u0,tspan,p=NullParameters();kwargs...)
+TwoPointBVProblem{isinplace}(f,bc,u0,tspan,p=NullParameters();kwargs...)
+BVProblem{isinplace}(f,bc,u0,tspan,p=NullParameters();kwargs...)
 ```
 
 or if we have an initial guess function `initialGuess(t)` for the given BVP,
 we can pass the initial guess to the problem constructors:
 
 ```julia
-TwoPointBVProblem{isinplace}(f,bc!,initialGuess,tspan,p=NullParameters();kwargs...)
-BVProblem{isinplace}(f,bc!,initialGuess,tspan,p=NullParameters();kwargs...)
+TwoPointBVProblem{isinplace}(f,bc,initialGuess,tspan,p=NullParameters();kwargs...)
+BVProblem{isinplace}(f,bc,initialGuess,tspan,p=NullParameters();kwargs...)
 ```
 
-For any BVP problem type, `bc!` is the inplace function:
+For any BVP problem type, `bc` must be inplace if `f` is inplace. Otherwise it must be
+out-of-place.
+
+If the bvp is a StandardBVProblem (also known as a Multi-Point BV Problem) it must define
+either of the following functions
 
 ```julia
 bc!(residual, u, p, t)
+residual = bc(u, p, t)
 ```
 
 where `residual` computed from the current `u`. `u` is an array of solution values
@@ -61,6 +71,16 @@ time points, and for shooting type methods `u=sol` the ODE solution.
 Note that all features of the `ODESolution` are present in this form.
 In both cases, the size of the residual matches the size of the initial condition.
 
+If the bvp is a TwoPointBVProblem it must define either of the following functions
+
+```julia
+bc!((resid_a, resid_b), (u_a, u_b), p)
+resid_a, resid_b = bc((u_a, u_b), p)
+```
+
+where `resid_a` and `resid_b` are the residuals at the two endpoints, `u_a` and `u_b` are
+the solution values at the two endpoints, and `p` are the parameters.
+
 Parameters are optional, and if not given, then a `NullParameters()` singleton
 will be used which will throw nice errors if you try to index non-existent
 parameters. Any extra keyword arguments are passed on to the solvers. For example,
@@ -88,10 +108,11 @@ struct BVProblem{uType, tType, isinplace, P, F, BF, PT, K} <:
     problem_type::PT
     kwargs::K
 
-    @add_kwonly function BVProblem{iip}(f::AbstractBVPFunction{iip}, bc, u0, tspan,
-        p = NullParameters(), problem_type = StandardBVProblem(); kwargs...) where {iip}
+    @add_kwonly function BVProblem{iip}(f::AbstractBVPFunction{iip, TP}, bc, u0, tspan,
+        p = NullParameters(); kwargs...) where {iip, TP}
         _tspan = promote_tspan(tspan)
         warn_paramtype(p)
+        problem_type = TP ? TwoPointBVProblem() : StandardBVProblem()
         return new{typeof(u0), typeof(_tspan), iip, typeof(p), typeof(f), typeof(bc),
             typeof(problem_type), typeof(kwargs)}(f, bc, u0, _tspan, p, problem_type,
             kwargs)
@@ -113,42 +134,34 @@ function BVProblem(f::AbstractBVPFunction, u0, tspan, p = NullParameters(); kwar
     return BVProblem{isinplace(f)}(f, f.bc, u0, tspan, p; kwargs...)
 end
 
-"""
-$(TYPEDEF)
-"""
-struct TwoPointBVPFunction{bF}
-    bc::bF
-end
-TwoPointBVPFunction(; bc = error("No argument `bc`")) = TwoPointBVPFunction(bc)
-function (f::TwoPointBVPFunction)(residuala, residualb, ua, ub, p)
-    return f.bc(residuala, residualb, ua, ub, p)
-end
-function (f::TwoPointBVPFunction)(residual::Tuple, u, p)
-    return f(residual[1], residual[2], u[1], u[end], p)
-end
+# This is mostly a fake stuct and isn't used anywhere
+# But we need it for function calls like TwoPointBVProblem{iip}(...) = ...
+struct TwoPointBVPFunction{iip} end
 
-"""
-$(TYPEDEF)
-"""
-struct TwoPointBVProblem end
+@inline TwoPointBVPFunction(args...; kwargs...) = BVPFunction(args...; kwargs..., twopoint=true)
+@inline function TwoPointBVPFunction{iip}(args...; kwargs...) where {iip}
+    return BVPFunction{iip}(args...; kwargs..., twopoint=true)
+end
 
-function TwoPointBVProblem(f, bc, u0, tspan, p = NullParameters(); kwargs...)
-    return TwoPointBVProblem(BVPFunction(f, bc), u0, tspan, p; kwargs...)
+function TwoPointBVProblem(f, bc, u0, tspan, p = NullParameters();
+    bcresid_prototype=nothing, kwargs...)
+    return TwoPointBVProblem(TwoPointBVPFunction(f, bc; bcresid_prototype), u0, tspan, p;
+        kwargs...)
 end
-function TwoPointBVProblem(f::AbstractBVPFunction, u0, tspan, p = NullParameters();
-    kwargs...)
-    iip = isinplace(f)
-    return BVProblem{iip}(f, f.bc, u0, tspan, p, TwoPointBVProblem(); kwargs...)
+function TwoPointBVProblem(f::AbstractBVPFunction{iip, twopoint}, u0, tspan,
+    p = NullParameters(); kwargs...) where {iip, twopoint}
+    @assert twopoint "`TwoPointBVProblem` can only be used with a `TwoPointBVPFunction`. Instead of using `BVPFunction`, use `TwoPointBVPFunction` or pass a kwarg `twopoint=true` during the construction of the `BVPFunction`."
+    return BVProblem{iip}(f, f.bc, u0, tspan, p; kwargs...)
 end
 
 # Allow previous timeseries solution
 function TwoPointBVProblem(f::AbstractODEFunction, bc, sol::T, tspan::Tuple,
-    p = NullParameters()) where {T <: AbstractTimeseriesSolution}
-    TwoPointBVProblem(f, bc, sol.u, tspan, p)
+    p = NullParameters(); kwargs...) where {T <: AbstractTimeseriesSolution}
+    return TwoPointBVProblem(f, bc, sol.u, tspan, p; kwargs...)
 end
 # Allow initial guess function for the initial guess
 function TwoPointBVProblem(f::AbstractODEFunction, bc, initialGuess, tspan::AbstractVector,
     p = NullParameters(); kwargs...)
     u0 = [initialGuess(i) for i in tspan]
-    TwoPointBVProblem(f, bc, u0, (tspan[1], tspan[end]), p)
+    return TwoPointBVProblem(f, bc, u0, (tspan[1], tspan[end]), p; kwargs...)
 end
diff --git a/src/scimlfunctions.jl b/src/scimlfunctions.jl
index 38edeb094..4c7c0caba 100644
--- a/src/scimlfunctions.jl
+++ b/src/scimlfunctions.jl
@@ -2126,8 +2126,7 @@ TruncatedStacktraces.@truncate_stacktrace OptimizationFunction 1 2
 """
 $(TYPEDEF)
 """
-abstract type AbstractBVPFunction{iip} <:
-              AbstractDiffEqFunction{iip} end
+abstract type AbstractBVPFunction{iip, twopoint} <: AbstractDiffEqFunction{iip} end
 
 @doc doc"""
     BVPFunction{iip,F,BF,TMM,Ta,Tt,TJ,BCTJ,JVP,VJP,JP,BCJP,SP,TW,TWt,TPJ,S,S2,S3,O,TCV,BCTCV} <: AbstractBVPFunction{iip,specialize}
@@ -2232,12 +2231,9 @@ For more details on this argument, see the ODEFunction documentation.
 
 The fields of the BVPFunction type directly match the names of the inputs.
 """
-struct BVPFunction{iip, specialize, F, BF, TMM, Ta, Tt, TJ, BCTJ, JVP, VJP, JP,
-    BCJP, SP, TW, TWt,
-    TPJ,
-    S, S2, S3, O, TCV, BCTCV,
-    SYS} <:
-       AbstractBVPFunction{iip}
+struct BVPFunction{iip, specialize, twopoint, F, BF, TMM, Ta, Tt, TJ, BCTJ, JVP, VJP,
+    JP, BCJP, BCRP, SP, TW, TWt, TPJ, S, S2, S3, O, TCV, BCTCV,
+    SYS} <: AbstractBVPFunction{iip, twopoint}
     f::F
     bc::BF
     mass_matrix::TMM
@@ -3817,7 +3813,7 @@ function OptimizationFunction{iip}(f, adtype::AbstractADType = NoAD();
         cons_expr, sys)
 end
 
-function BVPFunction{iip, specialize}(f, bc;
+function BVPFunction{iip, specialize, twopoint}(f, bc;
     mass_matrix = __has_mass_matrix(f) ? f.mass_matrix : I,
     analytic = __has_analytic(f) ? f.analytic : nothing,
     tgrad = __has_tgrad(f) ? f.tgrad : nothing,
@@ -3838,8 +3834,7 @@ function BVPFunction{iip, specialize}(f, bc;
     observed = __has_observed(f) ? f.observed : DEFAULT_OBSERVED,
     colorvec = __has_colorvec(f) ? f.colorvec : nothing,
     bccolorvec = __has_colorvec(bc) ? bc.colorvec : nothing,
-    sys = __has_sys(f) ? f.sys : nothing) where {iip, specialize}
-    # TODO: Add checks for bcresid_prototype for TwoPointBVPFunction
+    sys = __has_sys(f) ? f.sys : nothing) where {iip, specialize, twopoint}
     if mass_matrix === I && typeof(f) <: Tuple
         mass_matrix = ((I for i in 1:length(f))...,)
     end
@@ -3879,7 +3874,7 @@ function BVPFunction{iip, specialize}(f, bc;
         _bccolorvec = bccolorvec
     end
 
-    bciip = isinplace(bc, 4, "bc", iip)
+    bciip = !twopoint ? isinplace(bc, 4, "bc", iip) : isinplace(bc, 3, "bc", iip)
     jaciip = jac !== nothing ? isinplace(jac, 4, "jac", iip) : iip
     bcjaciip = bcjac !== nothing ? isinplace(bcjac, 4, "bcjac", bciip) : bciip
     tgradiip = tgrad !== nothing ? isinplace(tgrad, 4, "tgrad", iip) : iip
@@ -3899,6 +3894,12 @@ function BVPFunction{iip, specialize}(f, bc;
         throw(NonconformingFunctionsError(functions))
     end
 
+    if iip && twopoint
+        if bcresid_prototype === nothing || length(bcresid_prototype) != 2
+            error("bcresid_prototype must be a tuple / indexable collection of length 2 for a TwoPointBVPFunction")
+        end
+    end
+
     if any(bc_nonconforming)
         bc_nonconforming = findall(bc_nonconforming)
         functions = ["bcjac"][bc_nonconforming]
@@ -3906,41 +3907,21 @@ function BVPFunction{iip, specialize}(f, bc;
     end
 
     if specialize === NoSpecialize
-        BVPFunction{iip, specialize, Any, Any, Any, Any, Any,
+        BVPFunction{iip, specialize, twopoint, Any, Any, Any, Any, Any,
             Any, Any, Any, Any, Any, Any, Any, Any, Any, Any,
             Any, typeof(syms), typeof(indepsym), typeof(paramsyms),
             Any, typeof(_colorvec), typeof(_bccolorvec), Any}(f, bc, mass_matrix,
-            analytic,
-            tgrad,
-            jac, bcjac, jvp, vjp,
-            jac_prototype,
+            analytic, tgrad, jac, bcjac, jvp, vjp, jac_prototype,
             bcjac_prototype, bcresid_prototype,
-            sparsity, Wfact,
-            Wfact_t,
-            paramjac, syms,
-            indepsym, paramsyms,
-            observed,
+            sparsity, Wfact, Wfact_t, paramjac, syms, indepsym, paramsyms, observed,
             _colorvec, _bccolorvec, sys)
-    elseif specialize === false
-        BVPFunction{iip, FunctionWrapperSpecialize,
-            typeof(f), typeof(bc), typeof(mass_matrix), typeof(analytic), typeof(tgrad),
-            typeof(jac), typeof(jvp), typeof(vjp), typeof(jac_prototype),
-            typeof(sparsity), typeof(Wfact), typeof(Wfact_t), typeof(paramjac),
-            typeof(syms), typeof(indepsym), typeof(paramsyms), typeof(observed),
-            typeof(_colorvec), typeof(_bccolorvec),
-            typeof(sys)}(f, bc, mass_matrix, analytic, tgrad, jac, bcjac,
-            jvp, vjp, jac_prototype, bcjac_prototype, sparsity, Wfact,
-            Wfact_t, paramjac, syms, indepsym, paramsyms,
-            observed, _colorvec, _bccolorvec, sys)
     else
-        BVPFunction{iip, specialize, typeof(f), typeof(bc), typeof(mass_matrix),
-            typeof(analytic),
-            typeof(tgrad),
-            typeof(jac), typeof(bcjac), typeof(jvp), typeof(vjp), typeof(jac_prototype),
-            typeof(bcjac_prototype), typeof(bcresid_prototype),
-            typeof(sparsity), typeof(Wfact), typeof(Wfact_t),
-            typeof(paramjac), typeof(syms), typeof(indepsym), typeof(paramsyms),
-            typeof(observed),
+        BVPFunction{iip, specialize, twopoint, typeof(f), typeof(bc), typeof(mass_matrix),
+            typeof(analytic), typeof(tgrad), typeof(jac), typeof(bcjac), typeof(jvp),
+            typeof(vjp), typeof(jac_prototype),
+            typeof(bcjac_prototype), typeof(bcresid_prototype), typeof(sparsity),
+            typeof(Wfact), typeof(Wfact_t), typeof(paramjac), typeof(syms),
+            typeof(indepsym), typeof(paramsyms), typeof(observed),
             typeof(_colorvec), typeof(_bccolorvec), typeof(sys)}(f, bc, mass_matrix, analytic,
             tgrad, jac, bcjac, jvp, vjp,
             jac_prototype, bcjac_prototype, bcresid_prototype, sparsity,
@@ -3950,12 +3931,12 @@ function BVPFunction{iip, specialize}(f, bc;
     end
 end
 
-function BVPFunction{iip}(f, bc; kwargs...) where {iip}
-    BVPFunction{iip, FullSpecialize}(f, bc; kwargs...)
+function BVPFunction{iip}(f, bc; twopoint::Bool=false, kwargs...) where {iip}
+    BVPFunction{iip, FullSpecialize, twopoint}(f, bc; kwargs...)
 end
 BVPFunction{iip}(f::BVPFunction, bc; kwargs...) where {iip} = f
-function BVPFunction(f, bc; kwargs...)
-    BVPFunction{isinplace(f, 4), FullSpecialize}(f, bc; kwargs...)
+function BVPFunction(f, bc; twopoint::Bool=false, kwargs...)
+    BVPFunction{isinplace(f, 4), FullSpecialize, twopoint}(f, bc; kwargs...)
 end
 BVPFunction(f::BVPFunction; kwargs...) = f
 

From cb2ad9fc1fe124983eb9b98f002557ba3d25051c Mon Sep 17 00:00:00 2001
From: Avik Pal <avikpal@mit.edu>
Date: Tue, 12 Sep 2023 13:52:02 -0400
Subject: [PATCH 06/16] Construct ArrayPartition

---
 Project.toml                 |  2 +-
 src/problems/bvp_problems.jl |  2 +-
 src/scimlfunctions.jl        | 15 +++++----------
 3 files changed, 7 insertions(+), 12 deletions(-)

diff --git a/Project.toml b/Project.toml
index 951f6b7ca..7b8a2b09a 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,7 +1,7 @@
 name = "SciMLBase"
 uuid = "0bca4576-84f4-4d90-8ffe-ffa030f20462"
 authors = ["Chris Rackauckas <accounts@chrisrackauckas.com> and contributors"]
-version = "2.0.0"
+version = "1.96.0"  # Remember to make it 2.0.0
 
 [deps]
 ADTypes = "47edcb42-4c32-4615-8424-f2b9edc5f35b"
diff --git a/src/problems/bvp_problems.jl b/src/problems/bvp_problems.jl
index b46b0479c..b6fe87866 100644
--- a/src/problems/bvp_problems.jl
+++ b/src/problems/bvp_problems.jl
@@ -127,7 +127,7 @@ TruncatedStacktraces.@truncate_stacktrace BVProblem 3 1 2
 
 function BVProblem(f, bc, u0, tspan, p = NullParameters(); kwargs...)
     iip = isinplace(f, 4)
-    return BVProblem{iip}(BVPFunction{iip}(f, bc), u0, tspan, p; kwargs...)
+    return BVProblem{iip}(BVPFunction{iip}(f, bc), bc, u0, tspan, p; kwargs...)
 end
 
 function BVProblem(f::AbstractBVPFunction, u0, tspan, p = NullParameters(); kwargs...)
diff --git a/src/scimlfunctions.jl b/src/scimlfunctions.jl
index 4c7c0caba..f7f2178bb 100644
--- a/src/scimlfunctions.jl
+++ b/src/scimlfunctions.jl
@@ -1884,10 +1884,7 @@ For more details on this argument, see the ODEFunction documentation.
 The fields of the NonlinearFunction type directly match the names of the inputs.
 """
 struct NonlinearFunction{iip, specialize, F, TMM, Ta, Tt, TJ, JVP, VJP, JP, SP, TW, TWt,
-    TPJ,
-    S, S2, O, TCV,
-    SYS,
-} <: AbstractNonlinearFunction{iip}
+    TPJ, S, S2, O, TCV, SYS} <: AbstractNonlinearFunction{iip}
     f::F
     mass_matrix::TMM
     analytic::Ta
@@ -3646,10 +3643,7 @@ function NonlinearFunction{iip, specialize}(f;
                DEFAULT_OBSERVED_NO_TIME,
     colorvec = __has_colorvec(f) ? f.colorvec :
                nothing,
-    sys = __has_sys(f) ? f.sys : nothing) where {
-    iip,
-    specialize,
-}
+    sys = __has_sys(f) ? f.sys : nothing) where {iip, specialize}
     if mass_matrix === I && typeof(f) <: Tuple
         mass_matrix = ((I for i in 1:length(f))...,)
     end
@@ -3701,8 +3695,8 @@ function NonlinearFunction{iip, specialize}(f;
             typeof(Wfact_t), typeof(paramjac), typeof(syms),
             typeof(paramsyms),
             typeof(observed),
-            typeof(_colorvec), typeof(sys)}(f, mass_matrix, analytic, tgrad,
-            jac,
+            typeof(_colorvec), typeof(sys)}(f, mass_matrix,
+            analytic, tgrad, jac,
             jvp, vjp, jac_prototype, sparsity,
             Wfact,
             Wfact_t, paramjac, syms,
@@ -3898,6 +3892,7 @@ function BVPFunction{iip, specialize, twopoint}(f, bc;
         if bcresid_prototype === nothing || length(bcresid_prototype) != 2
             error("bcresid_prototype must be a tuple / indexable collection of length 2 for a TwoPointBVPFunction")
         end
+        bcresid_prototype = ArrayPartition(bcresid_prototype[1], bcresid_prototype[2])
     end
 
     if any(bc_nonconforming)

From 8b006588b62c420262a5ff212a88f9dcb36abec7 Mon Sep 17 00:00:00 2001
From: Avik Pal <avikpal@mit.edu>
Date: Tue, 12 Sep 2023 18:07:17 -0400
Subject: [PATCH 07/16] Proper resid prototype handling

---
 src/scimlfunctions.jl | 10 ++++++----
 1 file changed, 6 insertions(+), 4 deletions(-)

diff --git a/src/scimlfunctions.jl b/src/scimlfunctions.jl
index 644545948..9da24fb14 100644
--- a/src/scimlfunctions.jl
+++ b/src/scimlfunctions.jl
@@ -3892,11 +3892,13 @@ function BVPFunction{iip, specialize, twopoint}(f, bc;
         throw(NonconformingFunctionsError(functions))
     end
 
-    if iip && twopoint
-        if bcresid_prototype === nothing || length(bcresid_prototype) != 2
-            error("bcresid_prototype must be a tuple / indexable collection of length 2 for a TwoPointBVPFunction")
+    if twopoint
+        if iip && (bcresid_prototype === nothing || length(bcresid_prototype) != 2)
+            error("bcresid_prototype must be a tuple / indexable collection of length 2 for a inplace TwoPointBVPFunction")
+        end
+        if bcresid_prototype !== nothing && length(bcresid_prototype) == 2
+            bcresid_prototype = ArrayPartition(bcresid_prototype[1], bcresid_prototype[2])
         end
-        bcresid_prototype = ArrayPartition(bcresid_prototype[1], bcresid_prototype[2])
     end
 
     if any(bc_nonconforming)

From ac73376f443bc832071ee54ee257ea978ac5db22 Mon Sep 17 00:00:00 2001
From: IlianPihlajamaa <73794090+IlianPihlajamaa@users.noreply.github.com>
Date: Mon, 18 Sep 2023 21:14:40 +0200
Subject: [PATCH 08/16] remove Value type requirement from
 SampledIntegralProblem

---
 src/problems/basic_problems.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/problems/basic_problems.jl b/src/problems/basic_problems.jl
index e0385140d..1dc56ec9f 100644
--- a/src/problems/basic_problems.jl
+++ b/src/problems/basic_problems.jl
@@ -437,7 +437,7 @@ struct SampledIntegralProblem{Y, X, D, K} <: AbstractIntegralProblem{false}
         @assert dim <= ndims(y) "The integration dimension `dim` is larger than the number of dimensions of the integrand `y`"
         @assert length(x)==size(y, dim) "The integrand `y` must have the same length as the sampling points `x` along the integrated dimension."
         @assert axes(x, 1)==axes(y, dim) "The integrand `y` must obey the same indexing as the sampling points `x` along the integrated dimension."
-        new{typeof(y), typeof(x), Val{dim}, typeof(kwargs)}(y, x, Val(dim), kwargs)
+        new{typeof(y), typeof(x), typeof(dim)}, typeof(kwargs)}(y, x, dim, kwargs)
     end
 end
 

From 819ec920ac894040c0687de94bbcc97b419a98d4 Mon Sep 17 00:00:00 2001
From: Christopher Rackauckas <accounts@chrisrackauckas.com>
Date: Tue, 19 Sep 2023 03:08:44 -0400
Subject: [PATCH 09/16] Update src/problems/basic_problems.jl

---
 src/problems/basic_problems.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/problems/basic_problems.jl b/src/problems/basic_problems.jl
index 1dc56ec9f..efa1bf3a0 100644
--- a/src/problems/basic_problems.jl
+++ b/src/problems/basic_problems.jl
@@ -437,7 +437,7 @@ struct SampledIntegralProblem{Y, X, D, K} <: AbstractIntegralProblem{false}
         @assert dim <= ndims(y) "The integration dimension `dim` is larger than the number of dimensions of the integrand `y`"
         @assert length(x)==size(y, dim) "The integrand `y` must have the same length as the sampling points `x` along the integrated dimension."
         @assert axes(x, 1)==axes(y, dim) "The integrand `y` must obey the same indexing as the sampling points `x` along the integrated dimension."
-        new{typeof(y), typeof(x), typeof(dim)}, typeof(kwargs)}(y, x, dim, kwargs)
+        new{typeof(y), typeof(x), typeof(dim), typeof(kwargs)}(y, x, dim, kwargs)
     end
 end
 

From d70980ed563d985cc145923dc9249c11fc792190 Mon Sep 17 00:00:00 2001
From: Christopher Rackauckas <accounts@chrisrackauckas.com>
Date: Wed, 20 Sep 2023 16:59:10 -0400
Subject: [PATCH 10/16] Update basic_problems.jl

---
 src/problems/basic_problems.jl | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/src/problems/basic_problems.jl b/src/problems/basic_problems.jl
index efa1bf3a0..cd057eba5 100644
--- a/src/problems/basic_problems.jl
+++ b/src/problems/basic_problems.jl
@@ -426,10 +426,10 @@ SampledIntegralProblem(y::AbstractArray, x::AbstractVector; dim=ndims(y), kwargs
 
 The fields match the names of the constructor arguments.
 """
-struct SampledIntegralProblem{Y, X, D, K} <: AbstractIntegralProblem{false}
+struct SampledIntegralProblem{Y, X, K} <: AbstractIntegralProblem{false}
     y::Y
     x::X
-    dim::D
+    dim::Int
     kwargs::K
     @add_kwonly function SampledIntegralProblem(y::AbstractArray, x::AbstractVector;
         dim = ndims(y),
@@ -437,7 +437,7 @@ struct SampledIntegralProblem{Y, X, D, K} <: AbstractIntegralProblem{false}
         @assert dim <= ndims(y) "The integration dimension `dim` is larger than the number of dimensions of the integrand `y`"
         @assert length(x)==size(y, dim) "The integrand `y` must have the same length as the sampling points `x` along the integrated dimension."
         @assert axes(x, 1)==axes(y, dim) "The integrand `y` must obey the same indexing as the sampling points `x` along the integrated dimension."
-        new{typeof(y), typeof(x), typeof(dim), typeof(kwargs)}(y, x, dim, kwargs)
+        new{typeof(y), typeof(x), typeof(kwargs)}(y, x, dim, kwargs)
     end
 end
 

From eb9299515d9c8f50d08dc6406f2a929eb23d8a06 Mon Sep 17 00:00:00 2001
From: Christopher Rackauckas <accounts@chrisrackauckas.com>
Date: Wed, 20 Sep 2023 16:59:47 -0400
Subject: [PATCH 11/16] Update Project.toml

---
 Project.toml | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/Project.toml b/Project.toml
index fa982a687..95a458cb1 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,7 +1,7 @@
 name = "SciMLBase"
 uuid = "0bca4576-84f4-4d90-8ffe-ffa030f20462"
 authors = ["Chris Rackauckas <accounts@chrisrackauckas.com> and contributors"]
-version = "1.98.0"
+version = "1.98.1"
 
 [deps]
 ADTypes = "47edcb42-4c32-4615-8424-f2b9edc5f35b"

From ee4931ceb04730185b0a8c84d6c94bb2f1f4a7dc Mon Sep 17 00:00:00 2001
From: Avik Pal <avikpal@amdci2.julia.csail.mit.edu>
Date: Thu, 21 Sep 2023 13:28:45 -0400
Subject: [PATCH 12/16] Fix the remake issue with problem_type

---
 src/problems/bvp_problems.jl | 10 ++++++++--
 1 file changed, 8 insertions(+), 2 deletions(-)

diff --git a/src/problems/bvp_problems.jl b/src/problems/bvp_problems.jl
index b6fe87866..2e9f1feab 100644
--- a/src/problems/bvp_problems.jl
+++ b/src/problems/bvp_problems.jl
@@ -109,10 +109,16 @@ struct BVProblem{uType, tType, isinplace, P, F, BF, PT, K} <:
     kwargs::K
 
     @add_kwonly function BVProblem{iip}(f::AbstractBVPFunction{iip, TP}, bc, u0, tspan,
-        p = NullParameters(); kwargs...) where {iip, TP}
+        p = NullParameters(); problem_type=nothing, kwargs...) where {iip, TP}
         _tspan = promote_tspan(tspan)
         warn_paramtype(p)
-        problem_type = TP ? TwoPointBVProblem() : StandardBVProblem()
+        prob_type = TP ? TwoPointBVProblem() : StandardBVProblem()
+        # Needed to ensure that `problem_type` doesn't get passed in kwargs
+        if problem_type === nothing
+            problem_type = prob_type
+        else
+            @assert prob_type === problem_type "This indicates incorrect problem type specification! Users should never pass in `problem_type` kwarg, this exists exclusively for internal use."
+        end
         return new{typeof(u0), typeof(_tspan), iip, typeof(p), typeof(f), typeof(bc),
             typeof(problem_type), typeof(kwargs)}(f, bc, u0, _tspan, p, problem_type,
             kwargs)

From 7b2d7f74aa5d2b77d91b00aa997f4c497b321299 Mon Sep 17 00:00:00 2001
From: Avik Pal <avikpal@amdci2.julia.csail.mit.edu>
Date: Thu, 21 Sep 2023 13:44:02 -0400
Subject: [PATCH 13/16] Fix BVP tests

---
 Project.toml                             |  2 +-
 test/downstream/ensemble_bvp.jl          |  2 +-
 test/function_building_error_messages.jl | 49 +++++++++++++++---------
 3 files changed, 32 insertions(+), 21 deletions(-)

diff --git a/Project.toml b/Project.toml
index c63cead23..492af0df1 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,7 +1,7 @@
 name = "SciMLBase"
 uuid = "0bca4576-84f4-4d90-8ffe-ffa030f20462"
 authors = ["Chris Rackauckas <accounts@chrisrackauckas.com> and contributors"]
-version = "1.98.2"  # Remember to make it 2.0.0
+version = "2.0.0"
 
 [deps]
 ADTypes = "47edcb42-4c32-4615-8424-f2b9edc5f35b"
diff --git a/test/downstream/ensemble_bvp.jl b/test/downstream/ensemble_bvp.jl
index 4f3b2bc2a..ad08236c7 100644
--- a/test/downstream/ensemble_bvp.jl
+++ b/test/downstream/ensemble_bvp.jl
@@ -19,4 +19,4 @@ tspan = (0.0, pi / 2)
 p = [rand()]
 bvp = BVProblem(ode!, bc!, initial_guess, tspan, p)
 ensemble_prob = EnsembleProblem(bvp, prob_func = prob_func)
-sim = solve(ensemble_prob, GeneralMIRK4(), trajectories = 10, dt = 0.1)
+sim = solve(ensemble_prob, MIRK4(), trajectories = 10, dt = 0.1)
diff --git a/test/function_building_error_messages.jl b/test/function_building_error_messages.jl
index c424a6f84..279c80e73 100644
--- a/test/function_building_error_messages.jl
+++ b/test/function_building_error_messages.jl
@@ -463,7 +463,7 @@ IntegralProblem(intf, [0.0], [1.0], p)
 x = [1.0, 2.0]
 y = rand(2, 2)
 SampledIntegralProblem(y, x)
-SampledIntegralProblem(y, x; dim=2)
+SampledIntegralProblem(y, x; dim = 2)
 
 # Optimization
 
@@ -516,10 +516,7 @@ bcjac(u, p, t) = [1.0]
     bcoop,
     jac = bjac,
     bcjac = bcjac)
-@test_throws SciMLBase.NonconformingFunctionsError BVPFunction(bfiip,
-    bciip,
-    jac = bjac,
-    bcjac = bcjac)
+BVPFunction(bfiip, bciip, jac = bjac, bcjac = bcjac)
 @test_throws SciMLBase.NonconformingFunctionsError BVPFunction(bfoop,
     bciip,
     jac = bjac,
@@ -528,8 +525,14 @@ BVPFunction(bfoop, bcoop, jac = bjac)
 bjac(du, u, p, t) = [1.0]
 bcjac(du, u, p, t) = [1.0]
 BVPFunction(bfiip, bciip, jac = bjac, bcjac = bcjac)
-BVPFunction(bfoop, bciip, jac = bjac, bcjac = bcjac)
-BVPFunction(bfiip, bcoop, jac = bjac, bcjac = bcjac)
+@test_throws SciMLBase.NonconformingFunctionsError BVPFunction(bfoop,
+    bciip,
+    jac = bjac,
+    bcjac = bcjac)
+@test_throws SciMLBase.NonconformingFunctionsError BVPFunction(bfiip,
+    bcoop,
+    jac = bjac,
+    bcjac = bcjac)
 BVPFunction(bfoop, bcoop, jac = bjac, bcjac = bcjac)
 
 bWfact(u, t) = [1.0]
@@ -540,10 +543,10 @@ bWfact(u, p, t) = [1.0]
 @test_throws SciMLBase.TooFewArgumentsError BVPFunction(bfoop, bciip, Wfact = bWfact)
 bWfact(u, p, gamma, t) = [1.0]
 @test_throws SciMLBase.NonconformingFunctionsError BVPFunction(bfiip, bciip, Wfact = bWfact)
-BVPFunction(bfoop, bciip, Wfact = bWfact)
+@test_throws SciMLBase.NonconformingFunctionsError BVPFunction(bfoop, bciip, Wfact = bWfact)
 bWfact(du, u, p, gamma, t) = [1.0]
 BVPFunction(bfiip, bciip, Wfact = bWfact)
-BVPFunction(bfoop, bciip, Wfact = bWfact)
+@test_throws SciMLBase.NonconformingFunctionsError BVPFunction(bfoop, bciip, Wfact = bWfact)
 
 bWfact_t(u, t) = [1.0]
 @test_throws SciMLBase.TooFewArgumentsError BVPFunction(bfiip, bciip, Wfact_t = bWfact_t)
@@ -555,20 +558,24 @@ bWfact_t(u, p, gamma, t) = [1.0]
 @test_throws SciMLBase.NonconformingFunctionsError BVPFunction(bfiip,
     bciip,
     Wfact_t = bWfact_t)
-BVPFunction(bfoop, bciip, Wfact_t = bWfact_t)
+@test_throws SciMLBase.NonconformingFunctionsError BVPFunction(bfoop,
+    bciip,
+    Wfact_t = bWfact_t)
 bWfact_t(du, u, p, gamma, t) = [1.0]
 BVPFunction(bfiip, bciip, Wfact_t = bWfact_t)
-BVPFunction(bfoop, bciip, Wfact_t = bWfact_t)
+@test_throws SciMLBase.NonconformingFunctionsError BVPFunction(bfoop,
+    bciip,
+    Wfact_t = bWfact_t)
 
 btgrad(u, t) = [1.0]
 @test_throws SciMLBase.TooFewArgumentsError BVPFunction(bfiip, bciip, tgrad = btgrad)
 @test_throws SciMLBase.TooFewArgumentsError BVPFunction(bfoop, bciip, tgrad = btgrad)
 btgrad(u, p, t) = [1.0]
 @test_throws SciMLBase.NonconformingFunctionsError BVPFunction(bfiip, bciip, tgrad = btgrad)
-BVPFunction(bfoop, bciip, tgrad = btgrad)
+@test_throws SciMLBase.NonconformingFunctionsError BVPFunction(bfoop, bciip, tgrad = btgrad)
 btgrad(du, u, p, t) = [1.0]
 BVPFunction(bfiip, bciip, tgrad = btgrad)
-BVPFunction(bfoop, bciip, tgrad = btgrad)
+@test_throws SciMLBase.NonconformingFunctionsError BVPFunction(bfoop, bciip, tgrad = btgrad)
 
 bparamjac(u, t) = [1.0]
 @test_throws SciMLBase.TooFewArgumentsError BVPFunction(bfiip, bciip, paramjac = bparamjac)
@@ -577,27 +584,31 @@ bparamjac(u, p, t) = [1.0]
 @test_throws SciMLBase.NonconformingFunctionsError BVPFunction(bfiip,
     bciip,
     paramjac = bparamjac)
-BVPFunction(bfoop, bciip, paramjac = bparamjac)
+@test_throws SciMLBase.NonconformingFunctionsError BVPFunction(bfoop,
+    bciip,
+    paramjac = bparamjac)
 bparamjac(du, u, p, t) = [1.0]
 BVPFunction(bfiip, bciip, paramjac = bparamjac)
-BVPFunction(bfoop, bciip, paramjac = bparamjac)
+@test_throws SciMLBase.NonconformingFunctionsError BVPFunction(bfoop,
+    bciip,
+    paramjac = bparamjac)
 
 bjvp(u, p, t) = [1.0]
 @test_throws SciMLBase.TooFewArgumentsError BVPFunction(bfiip, bciip, jvp = bjvp)
 @test_throws SciMLBase.TooFewArgumentsError BVPFunction(bfoop, bciip, jvp = bjvp)
 bjvp(u, v, p, t) = [1.0]
 @test_throws SciMLBase.NonconformingFunctionsError BVPFunction(bfiip, bciip, jvp = bjvp)
-BVPFunction(bfoop, bciip, jvp = bjvp)
+@test_throws SciMLBase.NonconformingFunctionsError BVPFunction(bfoop, bciip, jvp = bjvp)
 bjvp(du, u, v, p, t) = [1.0]
 BVPFunction(bfiip, bciip, jvp = bjvp)
-BVPFunction(bfoop, bciip, jvp = bjvp)
+@test_throws SciMLBase.NonconformingFunctionsError BVPFunction(bfoop, bciip, jvp = bjvp)
 
 bvjp(u, p, t) = [1.0]
 @test_throws SciMLBase.TooFewArgumentsError BVPFunction(bfiip, bciip, vjp = bvjp)
 @test_throws SciMLBase.TooFewArgumentsError BVPFunction(bfoop, bciip, vjp = bvjp)
 bvjp(u, v, p, t) = [1.0]
 @test_throws SciMLBase.NonconformingFunctionsError BVPFunction(bfiip, bciip, vjp = bvjp)
-BVPFunction(bfoop, bciip, vjp = bvjp)
+@test_throws SciMLBase.NonconformingFunctionsError BVPFunction(bfoop, bciip, vjp = bvjp)
 bvjp(du, u, v, p, t) = [1.0]
 BVPFunction(bfiip, bciip, vjp = bvjp)
-BVPFunction(bfoop, bciip, vjp = bvjp)
+@test_throws SciMLBase.NonconformingFunctionsError BVPFunction(bfoop, bciip, vjp = bvjp)

From 8b173b96e74bab2b775981775b20d8c0b94bb138 Mon Sep 17 00:00:00 2001
From: Avik Pal <avikpal@mit.edu>
Date: Thu, 21 Sep 2023 13:58:08 -0400
Subject: [PATCH 14/16] Fix 1 last test

---
 test/function_building_error_messages.jl | 5 ++++-
 1 file changed, 4 insertions(+), 1 deletion(-)

diff --git a/test/function_building_error_messages.jl b/test/function_building_error_messages.jl
index 279c80e73..1db0fb468 100644
--- a/test/function_building_error_messages.jl
+++ b/test/function_building_error_messages.jl
@@ -516,7 +516,10 @@ bcjac(u, p, t) = [1.0]
     bcoop,
     jac = bjac,
     bcjac = bcjac)
-BVPFunction(bfiip, bciip, jac = bjac, bcjac = bcjac)
+@test_throws SciMLBase.NonconformingFunctionsError BVPFunction(bfiip,
+    bciip,
+    jac = bjac,
+    bcjac = bcjac)
 @test_throws SciMLBase.NonconformingFunctionsError BVPFunction(bfoop,
     bciip,
     jac = bjac,

From a64db934fbefe1e581d981f22dca9700c139b201 Mon Sep 17 00:00:00 2001
From: Avik Pal <avikpal@mit.edu>
Date: Thu, 21 Sep 2023 13:58:46 -0400
Subject: [PATCH 15/16] DONT MERGE TEST COMMIT

---
 Project.toml | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/Project.toml b/Project.toml
index 492af0df1..b15c4569b 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,7 +1,7 @@
 name = "SciMLBase"
 uuid = "0bca4576-84f4-4d90-8ffe-ffa030f20462"
 authors = ["Chris Rackauckas <accounts@chrisrackauckas.com> and contributors"]
-version = "2.0.0"
+version = "1.99.0"
 
 [deps]
 ADTypes = "47edcb42-4c32-4615-8424-f2b9edc5f35b"

From 593e7f05ee51aea1800e86aed0b25f54de5eb17f Mon Sep 17 00:00:00 2001
From: Avik Pal <avikpal@mit.edu>
Date: Thu, 21 Sep 2023 14:33:20 -0400
Subject: [PATCH 16/16] Bump Project version

---
 Project.toml | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/Project.toml b/Project.toml
index b15c4569b..492af0df1 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,7 +1,7 @@
 name = "SciMLBase"
 uuid = "0bca4576-84f4-4d90-8ffe-ffa030f20462"
 authors = ["Chris Rackauckas <accounts@chrisrackauckas.com> and contributors"]
-version = "1.99.0"
+version = "2.0.0"
 
 [deps]
 ADTypes = "47edcb42-4c32-4615-8424-f2b9edc5f35b"