diff --git a/lib/BVProblemLibrary/src/linear.jl b/lib/BVProblemLibrary/src/linear.jl
index 0947e39..f368cc9 100644
--- a/lib/BVProblemLibrary/src/linear.jl
+++ b/lib/BVProblemLibrary/src/linear.jl
@@ -20,6 +20,22 @@ end
 prob_bvp_linear_1_function = BVPFunction(
     prob_bvp_linear_1_f!, (prob_bvp_linear_1_bca!, prob_bvp_linear_1_bcb!),
     analytic = prob_bvp_linear_1_analytic, bcresid_prototype = (zeros(1), zeros(1)), twopoint = Val(true))
+
+#### Out Of Place (OOP) version
+function prob_bvp_linear_1_f_oop(u, p, t)
+    du1 = u[2]
+    du2 = 1 / p * u[1]
+    return [du1, du2]
+end
+function prob_bvp_linear_1_bca(u_a, p)
+    return u_a[1] - 1
+end
+function prob_bvp_linear_1_bcb(u_b, p)
+    return u_b[1]
+end
+prob_bvp_linear_1_function_oop = BVPFunction(
+    prob_bvp_linear_1_f_oop, @SVector [prob_bvp_linear_1_bca; prob_bvp_linear_1_bcb],
+        analytic = prob_bvp_linear_1_analytic, bcresid_prototype = @SVector [zeros(1), zeros(1)], twopoint = Val(true))
 prob_bvp_linear_1_tspan = (0.0, 1.0)
 @doc raw"""
     prob_bvp_linear_1
@@ -59,6 +75,11 @@ prob_bvp_linear_1 = BVProblem(prob_bvp_linear_1_function,
     prob_bvp_linear_1_tspan,
     λ)
 
+prob_bvp_linear_1_oop = BVProblem(prob_bvp_linear_1_function_oop,
+    [1.0, 0.0],
+    prob_bvp_linear_1_tspan,
+    λ)
+
 ################### linear_bvp2 ############################
 function prob_bvp_linear_2_analytic(u, λ, t)
     a = 1 / sqrt(λ)
@@ -78,6 +99,22 @@ end
 prob_bvp_linear_2_function = BVPFunction(
     prob_bvp_linear_2_f!, (prob_bvp_linear_2_bca!, prob_bvp_linear_2_bcb!),
     analytic = prob_bvp_linear_2_analytic, bcresid_prototype = (zeros(1), zeros(1)), twopoint = Val(true))
+
+### OOP version
+function prob_bvp_linear_2_f_oop(u, p, t)
+    du1 = u[2]
+    du2 = 1 / p * u[2]
+    return [du1, du2]
+end
+function prob_bvp_linear_2_bca(u_a, p)
+    return u_a[1] - 1
+end
+function prob_bvp_linear_2_bcb(u_b, p)
+    return u_b[1]
+end
+prob_bvp_linear_2_function_oop = BVPFunction(
+    prob_bvp_linear_2_f_oop, @SVector [prob_bvp_linear_2_bca; prob_bvp_linear_2_bcb],
+    analytic = prob_bvp_linear_2_analytic, bcresid_prototype = @SVector [zeros(1), zeros(1)], twopoint = Val(true))
 prob_bvp_linear_2_tspan = (0.0, 1.0)
 @doc raw"""
     prob_bvp_linear_2
@@ -116,6 +153,10 @@ prob_bvp_linear_2 = BVProblem(prob_bvp_linear_2_function,
     [1.0, 0.0],
     prob_bvp_linear_2_tspan,
     λ)
+prob_bvp_linear_2_oop = BVProblem(prob_bvp_linear_2_function_oop,
+    [1.0, 0.0],
+    prob_bvp_linear_2_tspan,
+    λ)
 
 ################### linear_bvp3 ############################
 function prob_bvp_linear_3_analytic(u, λ, t)
@@ -139,6 +180,23 @@ end
 prob_bvp_linear_3_function = BVPFunction(
     prob_bvp_linear_3_f!, (prob_bvp_linear_3_bca!, prob_bvp_linear_3_bcb!),
     bcresid_prototype = (zeros(1), zeros(1)), analytic = prob_bvp_linear_3_analytic, twopoint = Val(true))
+
+# OOP 
+function prob_bvp_linear_3_f_oop(u, p, t)
+    du1 = u[2]
+    du2 = 1 / p * prob_bvp_linear_3_f(t, u[1], u[2])
+    return [du1, du2]
+end
+function prob_bvp_linear_3_bca(u_a, p)
+    u_a[1] + 1
+end
+function prob_bvp_linear_3_bcb(u_b, p)
+    u_b[1] + 1
+end
+prob_bvp_linear_3_function_oop = BVPFunction(
+    prob_bvp_linear_3_f_oop, (prob_bvp_linear_3_bca, prob_bvp_linear_3_bcb),
+    bcresid_prototype = (zeros(1), zeros(1)), analytic = prob_bvp_linear_3_analytic, twopoint = Val(true))
+
 prob_bvp_linear_3_tspan = (-1.0, 1.0)
 @doc raw"""
     prob_bvp_linear_3
@@ -183,6 +241,10 @@ prob_bvp_linear_3 = BVProblem(prob_bvp_linear_3_function,
     [1.0, 0.0],
     prob_bvp_linear_3_tspan,
     λ)
+prob_bvp_linear_3_oop = BVProblem(prob_bvp_linear_3_function_oop,
+    [1.0, 0.0],
+    prob_bvp_linear_3_tspan,
+    λ)
 
 ################### linear_bvp4 ############################
 function prob_bvp_linear_4_analytic(u, λ, t)
@@ -203,6 +265,23 @@ end
 prob_bvp_linear_4_function = BVPFunction(
     prob_bvp_linear_4_f!, (prob_bvp_linear_4_bca!, prob_bvp_linear_4_bcb!),
     bcresid_prototype = (zeros(1), zeros(1)), analytic = prob_bvp_linear_4_analytic, twopoint = Val(true))
+
+#### OOP
+function prob_bvp_linear_4_f_oop(u, p, t)
+    du1 = u[2]
+    du2 = 1 / p * prob_bvp_linear_4_f(u[1], u[2], p)
+    return [du1, du2]
+end
+function prob_bvp_linear_4_bca(u_a, p)
+    u_a[1] - 1 - exp(-2)
+end
+function prob_bvp_linear_4_bcb(u_b, p)
+    u_b[1] - 1 - exp(-2 * (1 + p) / p)
+end
+prob_bvp_linear_4_function_oop = BVPFunction(
+    prob_bvp_linear_4_f_oop, (prob_bvp_linear_4_bca, prob_bvp_linear_4_bcb),
+    bcresid_prototype = (zeros(1), zeros(1)), analytic = prob_bvp_linear_4_analytic, twopoint = Val(true))
+
 prob_bvp_linear_4_tspan = (-1.0, 1.0)
 @doc raw"""
     prob_bvp_linear_4
@@ -247,6 +326,10 @@ prob_bvp_linear_4 = BVProblem(prob_bvp_linear_4_function,
     [1.0, 0.0],
     prob_bvp_linear_4_tspan,
     λ)
+prob_bvp_linear_4_oop = BVProblem(prob_bvp_linear_4_function_oop,
+    [1.0, 0.0],
+    prob_bvp_linear_4_tspan,
+    λ)
 
 ################### linear_bvp5 ############################
 function prob_bvp_linear_5_analytic(u, λ, t)
@@ -269,6 +352,23 @@ end
 prob_bvp_linear_5_function = BVPFunction(
     prob_bvp_linear_5_f!, (prob_bvp_linear_5_bca!, prob_bvp_linear_5_bcb!),
     bcresid_prototype = (zeros(1), zeros(1)), analytic = prob_bvp_linear_5_analytic, twopoint = Val(true))
+
+# OOP 
+function prob_bvp_linear_5_f_oop(u, p, t)
+    du1 = u[2]
+    du2 = 1 / p * prob_bvp_linear_5_f(t, u[1], u[2], p)
+    return [du1, du2]
+end
+function prob_bvp_linear_5_bca(u_a, p)
+    u_a[1] + 1
+end
+function prob_bvp_linear_5_bcb(u_b, p)
+    u_b[1] + 1
+end
+prob_bvp_linear_5_function_oop = BVPFunction(
+    prob_bvp_linear_5_f_oop, (prob_bvp_linear_5_bca, prob_bvp_linear_5_bcb),
+    bcresid_prototype = (zeros(1), zeros(1)), analytic = prob_bvp_linear_5_analytic, twopoint = Val(true))
+
 prob_bvp_linear_5_tspan = (-1.0, 1.0)
 @doc raw"""
     prob_bvp_linear_5
@@ -313,6 +413,10 @@ prob_bvp_linear_5 = BVProblem(prob_bvp_linear_5_function,
     [1.0, 0.0],
     prob_bvp_linear_5_tspan,
     λ)
+prob_bvp_linear_5_oop = BVProblem(prob_bvp_linear_5_function_oop,
+    [1.0, 0.0],
+    prob_bvp_linear_5_tspan,
+    λ)
 
 ################### linear_bvp6 ############################
 function prob_bvp_linear_6_analytic(u, λ, t)
@@ -334,6 +438,24 @@ end
 prob_bvp_linear_6_function = BVPFunction(
     prob_bvp_linear_6_f!, (prob_bvp_linear_6_bca!, prob_bvp_linear_6_bcb!),
     bcresid_prototype = (zeros(1), zeros(1)), analytic = prob_bvp_linear_6_analytic, twopoint = Val(true))
+
+
+# OOP 
+function prob_bvp_linear_6_f_oop(u, p, t)
+    du1 = u[2]
+    du2 = 1 / p * prob_bvp_linear_6_f(t, u[2], p)
+    return [du1, du2]
+end
+function prob_bvp_linear_6_bca(u_a, p)
+    u_a[1] + 2
+end
+function prob_bvp_linear_6_bcb(u_b, p)
+    u_b[1]
+end
+prob_bvp_linear_6_function_oop = BVPFunction(
+    prob_bvp_linear_6_f_oop, (prob_bvp_linear_6_bca, prob_bvp_linear_6_bcb),
+    bcresid_prototype = (zeros(1), zeros(1)), analytic = prob_bvp_linear_6_analytic, twopoint = Val(true))
+
 prob_bvp_linear_6_tspan = (-1.0, 1.0)
 @doc raw"""
     prob_bvp_linear_6
@@ -378,6 +500,10 @@ prob_bvp_linear_6 = BVProblem(prob_bvp_linear_6_function,
     [1.0, 0.0],
     prob_bvp_linear_6_tspan,
     λ)
+prob_bvp_linear_6_oop = BVProblem(prob_bvp_linear_6_function_oop,
+    [1.0, 0.0],
+    prob_bvp_linear_6_tspan,
+    λ)
 
 ################### linear_bvp7 ############################
 function prob_bvp_linear_7_analytic(u, λ, t)
@@ -406,6 +532,23 @@ end
 prob_bvp_linear_7_function = BVPFunction(
     prob_bvp_linear_7_f!, (prob_bvp_linear_7_bca!, prob_bvp_linear_7_bcb!),
     bcresid_prototype = (zeros(1), zeros(1)), analytic = prob_bvp_linear_7_analytic, twopoint = Val(true))
+
+### OOP 
+function prob_bvp_linear_7_f_oop(u, p, t)
+    du1 = u[2]
+    du2 = 1 / p * prob_bvp_linear_7_f(t, u[1], u[2], p)
+    return [du1, du2]
+end
+function prob_bvp_linear_7_bca(u_a, p)
+    u_a[1] + 1
+end
+function prob_bvp_linear_7_bcb(u_b, p)
+    u_b[1] - 1
+end
+prob_bvp_linear_7_function_oop = BVPFunction(
+    prob_bvp_linear_7_f_oop, (prob_bvp_linear_7_bca, prob_bvp_linear_7_bcb),
+    bcresid_prototype = (zeros(1), zeros(1)), analytic = prob_bvp_linear_7_analytic, twopoint = Val(true))
+
 prob_bvp_linear_7_tspan = (-1.0, 1.0)
 @doc raw"""
     prob_bvp_linear_7
@@ -450,6 +593,10 @@ prob_bvp_linear_7 = BVProblem(prob_bvp_linear_7_function,
     [1.0, 0.0],
     prob_bvp_linear_7_tspan,
     λ)
+prob_bvp_linear_7_oop = BVProblem(prob_bvp_linear_7_function_oop,
+    [1.0, 0.0],
+    prob_bvp_linear_7_tspan,
+    λ)
 
 ################### linear_bvp8 ############################
 function prob_bvp_linear_8_analytic(u, λ, t)
@@ -470,6 +617,24 @@ end
 prob_bvp_linear_8_function = BVPFunction(
     prob_bvp_linear_8_f!, (prob_bvp_linear_8_bca!, prob_bvp_linear_8_bcb!),
     bcresid_prototype = (zeros(1), zeros(1)), analytic = prob_bvp_linear_8_analytic, twopoint = Val(true))
+
+# OOP
+function prob_bvp_linear_8_f_oop(u, p, t)
+    du1 = u[2]
+    du2 = -1 / p * u[2]
+    return [du1, du2]
+end
+function prob_bvp_linear_8_bca(u_a, p)
+    u_a[1] - 1
+end
+function prob_bvp_linear_8_bcb(u_b, p)
+    u_b[1] - 2
+end
+prob_bvp_linear_8_function_oop = BVPFunction(
+    prob_bvp_linear_8_f_oop, (prob_bvp_linear_8_bca, prob_bvp_linear_8_bcb),
+    bcresid_prototype = (zeros(1), zeros(1)), analytic = prob_bvp_linear_8_analytic, twopoint = Val(true))
+
+
 prob_bvp_linear_8_tspan = (0.0, 1.0)
 @doc raw"""
     prob_bvp_linear_8
@@ -508,6 +673,10 @@ prob_bvp_linear_8 = BVProblem(prob_bvp_linear_8_function,
     [1.0, 0.0],
     prob_bvp_linear_8_tspan,
     λ)
+prob_bvp_linear_8_oop = BVProblem(prob_bvp_linear_8_function_oop,
+    [1.0, 0.0],
+    prob_bvp_linear_8_tspan,
+    λ)
 
 ################### linear_bvp9 ############################
 function prob_bvp_linear_9_analytic(u, λ, t)
@@ -528,6 +697,23 @@ end
 prob_bvp_linear_9_function = BVPFunction(
     prob_bvp_linear_9_f!, (prob_bvp_linear_9_bca!, prob_bvp_linear_9_bcb!),
     bcresid_prototype = (zeros(1), zeros(1)), analytic = prob_bvp_linear_9_analytic, twopoint = Val(true))
+
+# OOP
+function prob_bvp_linear_9_f_oop(u, p, t)
+    du1 = u[2]
+    du2 = 1 / (p + t^2) * prob_bvp_linear_9_f(t, u[1], u[2], p)
+    return [du1, du2]
+end
+function prob_bvp_linear_9_bca(u_a, p)
+    u_a[1] - 1 / (1 + p)
+end
+function prob_bvp_linear_9_bcb(u_b, p)
+    u_b[1] - 1 / (1 + p)
+end
+prob_bvp_linear_9_function_oop = BVPFunction(
+    prob_bvp_linear_9_f_oop, (prob_bvp_linear_9_bca, prob_bvp_linear_9_bcb),
+    bcresid_prototype = (zeros(1), zeros(1)), analytic = prob_bvp_linear_9_analytic, twopoint = Val(true))
+
 prob_bvp_linear_9_tspan = (-1.0, 1.0)
 @doc raw"""
     prob_bvp_linear_9
@@ -572,6 +758,10 @@ prob_bvp_linear_9 = BVProblem(prob_bvp_linear_9_function,
     [1.0, 0.0],
     prob_bvp_linear_9_tspan,
     λ)
+prob_bvp_linear_9_oop = BVProblem(prob_bvp_linear_9_function_oop,
+    [1.0, 0.0],
+    prob_bvp_linear_9_tspan,
+    λ)
 
 ################### linear_bvp10 ############################
 function prob_bvp_linear_10_analytic(u, λ, t)
@@ -582,7 +772,7 @@ end
 prob_bvp_linear_10_f(t, u2, λ) = -t * u2
 function prob_bvp_linear_10_f!(du, u, p, t)
     du[1] = u[2]
-    du[2] = 1 / p * prob_bvp_linear_10_f(t, u[2], p)
+    du[2] = 1 / p * prob_bvp_linear_10f(t, u[2], p)
 end
 function prob_bvp_linear_10_bca!(res_a, u_a, p)
     res_a[1] = u_a[1]
@@ -593,6 +783,23 @@ end
 prob_bvp_linear_10_function = BVPFunction(
     prob_bvp_linear_10_f!, (prob_bvp_linear_10_bca!, prob_bvp_linear_10_bcb!),
     bcresid_prototype = (zeros(1), zeros(1)), analytic = prob_bvp_linear_10_analytic, twopoint = Val(true))
+
+# OOP
+function prob_bvp_linear_10_f_oop(u, p, t)
+    du1 = u[2]
+    du2 = 1 / p * prob_bvp_linear_10_f(t, u[2], p)
+    return [du1, du2]
+end
+function prob_bvp_linear_10_bca(u_a, p)
+    u_a[1]
+end
+function prob_bvp_linear_10_bcb(u_b, p)
+    u_b[1] - 2
+end
+prob_bvp_linear_10_function_oop = BVPFunction(
+    prob_bvp_linear_10_f_oop, (prob_bvp_linear_10_bca, prob_bvp_linear_10_bcb),
+    bcresid_prototype = (zeros(1), zeros(1)), analytic = prob_bvp_linear_10_analytic, twopoint = Val(true))
+
 prob_bvp_linear_10_tspan = (-1.0, 1.0)
 @doc raw"""
     prob_bvp_linear_10
@@ -637,6 +844,10 @@ prob_bvp_linear_10 = BVProblem(prob_bvp_linear_10_function,
     [1.0, 0.0],
     prob_bvp_linear_10_tspan,
     λ)
+prob_bvp_linear_10_oop = BVProblem(prob_bvp_linear_10_function_oop,
+    [1.0, 0.0],
+    prob_bvp_linear_10_tspan,
+    λ)
 
 ################### linear_bvp11 ############################
 function prob_bvp_linear_11_analytic(u, λ, t)
@@ -657,6 +868,23 @@ end
 prob_bvp_linear_11_function = BVPFunction(
     prob_bvp_linear_11_f!, (prob_bvp_linear_11_bca!, prob_bvp_linear_11_bcb!),
     bcresid_prototype = (zeros(1), zeros(1)), analytic = prob_bvp_linear_11_analytic, twopoint = Val(true))
+
+# OOP 
+function prob_bvp_linear_11_f_oop(u, p, t)
+    du1 = u[2]
+    du2 = 1 / p * prob_bvp_linear_11_f(t, u[1], p)
+    return [du1, du2]
+end
+function prob_bvp_linear_11_bca(u_a, p)
+    u_a[1] + 1
+end
+function prob_bvp_linear_11_bcb(u_b, p)
+    u_b[1] + 1
+end
+prob_bvp_linear_11_function_oop = BVPFunction(
+    prob_bvp_linear_11_f_oop, (prob_bvp_linear_11_bca, prob_bvp_linear_11_bcb),
+    bcresid_prototype = (zeros(1), zeros(1)), analytic = prob_bvp_linear_11_analytic, twopoint = Val(true))
+
 prob_bvp_linear_11_tspan = (-1.0, 1.0)
 @doc raw"""
     prob_bvp_linear_11
@@ -701,6 +929,10 @@ prob_bvp_linear_11 = BVProblem(prob_bvp_linear_11_function,
     [1.0, 0.0],
     prob_bvp_linear_11_tspan,
     λ)
+prob_bvp_linear_11_oop = BVProblem(prob_bvp_linear_11_function_oop,
+    [1.0, 0.0],
+    prob_bvp_linear_11_tspan,
+    λ)
 
 ################### linear_bvp12 ############################
 function prob_bvp_linear_12_analytic(u, λ, t)
@@ -712,7 +944,7 @@ end
 prob_bvp_linear_12_f(t, u1, λ) = u1 - (1 + λ * π^2) * cos(π * t)
 function prob_bvp_linear_12_f!(du, u, p, t)
     du[1] = u[2]
-    du[2] = 1 / p * prob_bvp_linear_12_f(t, u[1], p)
+    du[2] = 1 / p * prob_bvp_linear_12f(t, u[1], p)
 end
 function prob_bvp_linear_12_bca!(res_a, u_a, p)
     res_a[1] = u_a[1] + 1
@@ -723,6 +955,23 @@ end
 prob_bvp_linear_12_function = BVPFunction(
     prob_bvp_linear_12_f!, (prob_bvp_linear_12_bca!, prob_bvp_linear_12_bcb!),
     bcresid_prototype = (zeros(1), zeros(1)), analytic = prob_bvp_linear_12_analytic, twopoint = Val(true))
+
+# OOP 
+function prob_bvp_linear_12_f_oop(u, p, t)
+    du1 = u[2]
+    du2 = 1 / p * prob_bvp_linear_12_f(t, u[1], p)
+    return [du1, du2]
+end
+function prob_bvp_linear_12_bca(u_a, p)
+    u_a[1] + 1
+end
+function prob_bvp_linear_12_bcb(u_b, p)
+    u_b[1]
+end
+prob_bvp_linear_12_function_oop = BVPFunction(
+    prob_bvp_linear_12_f_oop, (prob_bvp_linear_12_bca, prob_bvp_linear_12_bcb),
+    bcresid_prototype = (zeros(1), zeros(1)), analytic = prob_bvp_linear_12_analytic, twopoint = Val(true))
+
 prob_bvp_linear_12_tspan = (-1.0, 1.0)
 @doc raw"""
     prob_bvp_linear_12
@@ -767,6 +1016,10 @@ prob_bvp_linear_12 = BVProblem(prob_bvp_linear_12_function,
     [1.0, 0.0],
     prob_bvp_linear_12_tspan,
     λ)
+prob_bvp_linear_12_oop = BVProblem(prob_bvp_linear_12_function_oop,
+    [1.0, 0.0],
+    prob_bvp_linear_12_tspan,
+    λ)
 
 ################### linear_bvp13 ############################
 function prob_bvp_linear_13_analytic(u, λ, t)
@@ -788,6 +1041,24 @@ end
 prob_bvp_linear_13_function = BVPFunction(
     prob_bvp_linear_13_f!, (prob_bvp_linear_13_bca!, prob_bvp_linear_13_bcb!),
     bcresid_prototype = (zeros(1), zeros(1)), analytic = prob_bvp_linear_13_analytic, twopoint = Val(true))
+
+# OOP 
+function prob_bvp_linear_13_f_oop(u, p, t)
+    du1 = u[2]
+    du2 = 1 / p * prob_bvp_linear_13_f(t, u[1], p)
+    return [du1, du2]
+end
+function prob_bvp_linear_13_bca(u_a, p)
+    u_a[1]
+end
+function prob_bvp_linear_13_bcb(u_b, p)
+    u_b[1] + 1 - exp(-2 / sqrt(p))
+end
+prob_bvp_linear_13_function_oop = BVPFunction(
+    prob_bvp_linear_13_f_oop, (prob_bvp_linear_13_bca, prob_bvp_linear_13_bcb),
+    bcresid_prototype = (zeros(1), zeros(1)), analytic = prob_bvp_linear_13_analytic, twopoint = Val(true))
+
+
 prob_bvp_linear_13_tspan = (-1.0, 1.0)
 @doc raw"""
     prob_bvp_linear_13
@@ -832,6 +1103,10 @@ prob_bvp_linear_13 = BVProblem(prob_bvp_linear_13_function,
     [1.0, 0.0],
     prob_bvp_linear_13_tspan,
     λ)
+prob_bvp_linear_13_oop = BVProblem(prob_bvp_linear_13_function_oop,
+    [1.0, 0.0],
+    prob_bvp_linear_13_tspan,
+    λ)
 
 ################### linear_bvp14 ############################
 function prob_bvp_linear_14_analytic(u, λ, t)
@@ -853,6 +1128,23 @@ end
 prob_bvp_linear_14_function = BVPFunction(
     prob_bvp_linear_14_f!, (prob_bvp_linear_14_bca!, prob_bvp_linear_14_bcb!),
     bcresid_prototype = (zeros(1), zeros(1)), analytic = prob_bvp_linear_14_analytic, twopoint = Val(true))
+
+# OOP 
+function prob_bvp_linear_14_f_oop(u, p, t)
+    du1 = u[2]
+    du2 = 1 / p * prob_bvp_linear_14_f(t, u[1], p)
+    return [du1, du2]
+end
+function prob_bvp_linear_14_bca(u_a, p)
+    u_a[1] - exp(-2 / sqrt(p))
+end
+function prob_bvp_linear_14_bcb(u_b, p)
+    u_b[1] - exp(-2 / sqrt(p))
+end
+prob_bvp_linear_14_function_oop = BVPFunction(
+    prob_bvp_linear_14_f_oop, (prob_bvp_linear_14_bca, prob_bvp_linear_14_bcb),
+    bcresid_prototype = (zeros(1), zeros(1)), analytic = prob_bvp_linear_14_analytic, twopoint = Val(true))
+
 prob_bvp_linear_14_tspan = (-1.0, 1.0)
 @doc raw"""
     prob_bvp_linear_14
@@ -897,6 +1189,10 @@ prob_bvp_linear_14 = BVProblem(prob_bvp_linear_14_function,
     [1.0, 0.0],
     prob_bvp_linear_14_tspan,
     λ)
+prob_bvp_linear_14 = BVProblem(prob_bvp_linear_14_function_oop,
+    [1.0, 0.0],
+    prob_bvp_linear_14_tspan,
+    λ)
 
 ################### linear_bvp15 ############################
 # No analytical solution
@@ -914,6 +1210,23 @@ end
 prob_bvp_linear_15_function = BVPFunction(
     prob_bvp_linear_15_f!, (prob_bvp_linear_15_bca!, prob_bvp_linear_15_bcb!),
     bcresid_prototype = (zeros(1), zeros(1)), twopoint = Val(true))
+
+# OOP
+function prob_bvp_linear_15_f_oop(u, p, t)
+    du1 = u[2]
+    du2 = 1 / p * prob_bvp_linear_15_f(t, u[1], p)
+    return [du1, du2]
+end
+function prob_bvp_linear_15_bca(u_a, p)
+    u_a[1] - 1
+end
+function prob_bvp_linear_15_bcb!(u_b, p)
+    u_b[1] - 1
+end
+prob_bvp_linear_15_function_oop = BVPFunction(
+    prob_bvp_linear_15_f_oop, (prob_bvp_linear_15_bca, prob_bvp_linear_15_bcb),
+    bcresid_prototype = (zeros(1), zeros(1)), twopoint = Val(true))
+
 prob_bvp_linear_15_tspan = (-1.0, 1.0)
 @doc raw"""
     prob_bvp_linear_15
@@ -951,6 +1264,10 @@ prob_bvp_linear_15 = BVProblem(prob_bvp_linear_15_function,
     [1.0, 0.0],
     prob_bvp_linear_15_tspan,
     λ)
+prob_bvp_linear_15_oop = BVProblem(prob_bvp_linear_15_function_oop,
+    [1.0, 0.0],
+    prob_bvp_linear_15_tspan,
+    λ)
 
 ################### linear_bvp16 ############################
 function prob_bvp_linear_16_analytic(u, λ, t)
@@ -971,6 +1288,24 @@ end
 prob_bvp_linear_16_function = BVPFunction(
     prob_bvp_linear_16_f!, (prob_bvp_linear_16_bca!, prob_bvp_linear_16_bcb!),
     bcresid_prototype = (zeros(1), zeros(1)), analytic = prob_bvp_linear_16_analytic, twopoint = Val(true))
+
+# OOP
+function prob_bvp_linear_16_f_oop(du, u, p, t)
+    du1 = u[2]
+    du2 = -1 / p^2 * π^2 * u[1] / 4
+    return [du1, du2]
+end
+function prob_bvp_linear_16_bca(u_a, p)
+    u_a[1]
+end
+function prob_bvp_linear_16_bcb(u_b, p)
+    u_b[1] - sin(π / (2 * p))
+end
+prob_bvp_linear_16_function_oop = BVPFunction(
+    prob_bvp_linear_16_f_oop, (prob_bvp_linear_16_bca, prob_bvp_linear_16_bcb),
+    bcresid_prototype = (zeros(1), zeros(1)), analytic = prob_bvp_linear_16_analytic, twopoint = Val(true))
+
+
 prob_bvp_linear_16_tspan = (0.0, 1.0)
 @doc raw"""
     prob_bvp_linear_16
@@ -1015,6 +1350,10 @@ prob_bvp_linear_16 = BVProblem(prob_bvp_linear_16_function,
     [1.0, 0.0],
     prob_bvp_linear_16_tspan,
     λ)
+prob_bvp_linear_16_oop = BVProblem(prob_bvp_linear_16_function_oop,
+    [1.0, 0.0],
+    prob_bvp_linear_16_tspan,
+    λ)
 
 ################### linear_bvp17 ############################
 function prob_bvp_linear_17_analytic(u, λ, t)
@@ -1034,6 +1373,23 @@ end
 prob_bvp_linear_17_function = BVPFunction(
     prob_bvp_linear_17_f!, (prob_bvp_linear_17_bca!, prob_bvp_linear_17_bcb!),
     bcresid_prototype = (zeros(1), zeros(1)), analytic = prob_bvp_linear_17_analytic, twopoint = Val(true))
+
+# OOP
+function prob_bvp_linear_17_f_oop(u, p, t)
+    du1 = u[2]
+    du2 = -3 * p * u[1] / (p + t^2)^2
+    return [du1, du2]
+end
+function prob_bvp_linear_17_bca(u_a, p)
+    u_a[1] + 0.1 / sqrt(p + 0.01)
+end
+function prob_bvp_linear_17_bcb(u_b, p)
+    u_b[1] - 0.1 / sqrt(p + 0.01)
+end
+prob_bvp_linear_17_function_oop = BVPFunction(
+    prob_bvp_linear_17_f_oop, (prob_bvp_linear_17_bca, prob_bvp_linear_17_bcb),
+    bcresid_prototype = (zeros(1), zeros(1)), analytic = prob_bvp_linear_17_analytic, twopoint = Val(true))
+
 prob_bvp_linear_17_tspan = (-0.1, 0.1)
 @doc raw"""
     prob_bvp_linear_17
@@ -1078,6 +1434,10 @@ prob_bvp_linear_17 = BVProblem(prob_bvp_linear_17_function,
     [1.0, 0.0],
     prob_bvp_linear_17_tspan,
     λ)
+prob_bvp_linear_17_oop = BVProblem(prob_bvp_linear_17_function_oop,
+    [1.0, 0.0],
+    prob_bvp_linear_17_tspan,
+    λ)
 
 ################### linear_bvp18 ############################
 function prob_bvp_linear_18_analytic(u, λ, t)
@@ -1098,6 +1458,23 @@ end
 prob_bvp_linear_18_function = BVPFunction(
     prob_bvp_linear_18_f!, (prob_bvp_linear_18_bca!, prob_bvp_linear_18_bcb!),
     bcresid_prototype = (zeros(1), zeros(1)), analytic = prob_bvp_linear_18_analytic, twopoint = Val(true))
+
+# OOP 
+function prob_bvp_linear_18_f_oop(u, p, t)
+    du1 = u[2]
+    du2 = -1 / p * u[2]
+    return [du1, du2]
+end
+function prob_bvp_linear_18_bca(u_a, p)
+    u_a[1] - 1
+end
+function prob_bvp_linear_18_bcb(u_b, p)
+    u_b[1] - exp(-1 / p)
+end
+prob_bvp_linear_18_function_oop = BVPFunction(
+    prob_bvp_linear_18_f_oop, (prob_bvp_linear_18_bca, prob_bvp_linear_18_bcb),
+    bcresid_prototype = (zeros(1), zeros(1)), analytic = prob_bvp_linear_18_analytic, twopoint = Val(true))
+
 prob_bvp_linear_18_tspan = (0, 1)
 @doc raw"""
     prob_bvp_linear_18
@@ -1142,3 +1519,7 @@ prob_bvp_linear_18 = BVProblem(prob_bvp_linear_18_function,
     [1.0, 0.0],
     prob_bvp_linear_18_tspan,
     λ)
+prob_bvp_linear_18_oop = BVProblem(prob_bvp_linear_18_function_oop,
+    [1.0, 0.0],
+    prob_bvp_linear_18_tspan,
+    λ)