diff --git a/Project.toml b/Project.toml
index 44c5fe4..916cab5 100644
--- a/Project.toml
+++ b/Project.toml
@@ -5,6 +5,7 @@ version = "2.6.1"
 
 [deps]
 ArgCheck = "dce04be8-c92d-5529-be00-80e4d2c0e197"
+BandedMatrices = "aae01518-5342-5314-be14-df237901396f"
 DifferentiationInterface = "a0c0ee7d-e4b9-4e03-894e-1c5f64a51d63"
 ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
@@ -22,6 +23,7 @@ ToeplitzMatrices = "c751599d-da0a-543b-9d20-d0a503d91d24"
 
 [compat]
 ArgCheck = "2"
+BandedMatrices = "1"
 DifferentiationInterface = "=0.5.5"
 ForwardDiff = "0.10"
 LinearAlgebra = "1"
diff --git a/src/Fronts.jl b/src/Fronts.jl
index 3b2e99e..04c40b4 100644
--- a/src/Fronts.jl
+++ b/src/Fronts.jl
@@ -13,6 +13,7 @@ using DifferentiationInterface: value_and_derivative,
                                 value_derivative_and_second_derivative, AutoForwardDiff
 using ArgCheck: @argcheck
 using StaticArrays: @SVector, @SMatrix
+using BandedMatrices: BandedMatrix
 using RecursiveArrayTools: ArrayPartition
 using PCHIPInterpolation: Interpolator, integrate
 import NumericalIntegration
diff --git a/src/finite.jl b/src/finite.jl
index 8f73491..6b32dd6 100644
--- a/src/finite.jl
+++ b/src/finite.jl
@@ -349,7 +349,8 @@ function solve(prob::FiniteReservoirProblem{<:DiffusionEquation{1}},
         end
     end
 
-    f! = ODEFunction(f!)
+    f! = ODEFunction(f!,
+        jac_prototype = BandedMatrix{eltype(u)}(undef, (length(u), length(u)), (1, 2)))
 
     steady_state = DiscreteCallback(
         (u, t, integrator) -> all(u[(begin + 1):end] .≈ u[end]),