I have tried in place format, but I get similar errors. Below is an example solving the diffusion equation in 1 dimension: \frac{du}{dt} = \frac{d^2u}{dz^2}. Replacing spatial derivatives with center differences results in ODEs with a tridiagonal jacobian. In this case, my specified jacobian is not necessary because Sundials can calculate the jacobian with finite differences - this is just for demonstrations sake. How should i format the band-storage jacobian `J`

, and which options should be selected, so that Sundials can use it.

```
using OrdinaryDiffEq
using Sundials
using ForwardDiff
using BandedMatrices
struct parameters{T<:AbstractFloat}
v_dep::T # lower boundary condition
PhiT::T # upper flux boundary condition
end
function rhs(dudt,u,p::parameters,t) # right-hand-side of ODEs
n = length(u)
h = 1.0/n # spatial grid size
for i in 2:(n-1)
dudt[i] = (1/h^2.0)*(u[i-1] - 2.0*u[i] + u[i+1])
end
# lower boundary
PhiG = -p.v_dep*u[1] # flux into the bottom layer
dudt[1] = (1.0/h^2.0)*(-u[1] + u[2]) + PhiG/h
# upper boundary
dudt[end] = (1.0/h^2.0)*(u[end-1] - u[end]) + p.PhiT/h
nothing
end;
function jacobian(J,u,p::parameters,t)
dudt = Array{Float64,1}(undef,length(u))
rhs1 = (dudt,uu) -> rhs(dudt,uu,p,0.0)
ForwardDiff.jacobian!(J,rhs1,dudt,u0)
nothing
end
n = 1000 # number of ODEs.
u0 = zeros(Float64,n) # initial condition
params = parameters(10.0,1.0) # boundary conditions
tspan = [0.0,10.0]
# Solve with NO SPECIFIED JACOBIAN
prob1 = ODEProblem(rhs,u0,tspan,params)
sol1 = solve(prob1,CVODE_BDF(linear_solver=:Band,jac_upper=1,jac_lower=1)); # this works
# Solving WITH specified jacobian
J = BandedMatrix(zeros(Float64,(n,n)), (1,1))
jacobian(J,u0,params,0.0)
f = ODEFunction(rhs, jac=jacobian,jac_prototype=J)
prob2 = ODEProblem(f,u0,tspan,params)
sol2 = solve(prob2,Rodas5()); # this works
sol3 = solve(prob2,CVODE_BDF(linear_solver=:Band,jac_upper=1,jac_lower=1)) # Method error. CVODE can't accept BandedMatrix
```