I am trying to solve the following problem

in a one dimensional x\in[0, 1] domain, where p_1 and p_2 are parameters. The boundary conditions are

I was successful in implementing the forward solver using the `DAEProblem`

class.

```
using DiffEqOperators, OrdinaryDiffEq, Plots, DifferentialEquations, Zygote
nknots = 10
h = 1.0 / (nknots + 1)
knots = range(h, step=h, length=nknots)
ord_deriv = 2
ord_approx = 2
const Δ = CenteredDifference(ord_deriv, ord_approx, h, nknots)
τ = 0.1
ξ = 0.5
t0 = 0.0
tlimit = 1.0 / ξ
t1 = 2.0 * tlimit
vlimit = 1.0
p = [1.0, 2.0]
function heat!(res, du, u, p, t)
rbc_1 = RobinBC((1.0, 0.0, t), (0.0, 1.0, 0.0), h, 1);
rbc_2 = RobinBC((0.0, 1.0, 0.0), (1.0, 0.0, 0.0), h, 1);
res[:, 1] = du[:, 1] - du[:, 2] - 1.0 * Δ * rbc_1 * u[:, 1]
res[:, 2] = τ * p[1] * p[2] * (du[:, 1] - du[:, 2]) + 1.0 * Δ * rbc_2 * u[:, 2]
end
u0 = zeros(nknots, 2)
diff_var = trues(nknots, 2)
du0 = zeros(nknots, 2)
prob = DAEProblem{true}(heat!, du0, u0, (t0, t1)) # , differential_vars=diff_var)
sol = solve(prob, DABDF2(), abstol=1e-3, reltol=1e-3)
```

Now I’d like to get gradients with respect to p_1 and p_2, with the following lines:

```
function sum_of_solution(u0, p)
_prob = remake(prob, u0=u0, p=p)
sum(solve(_prob, DABDF2(), rtol=1e-6, atol=1e-6, saveat=0.1, sensealg=ReverseDiffAdjoint()))
end
du01, dp1 = Zygote.gradient(sum_of_solution, u0, p)
```

but I am getting the error:

```
ERROR: LoadError: MethodError: no method matching DiffEqOperators.BoundaryPaddedVector(::ReverseDiff.TrackedReal{Float64, Float64, Nothing}, ::ReverseDiff.TrackedReal{Float64, Float64, Nothing}, ::Vector{ReverseDiff.TrackedReal{Float64, Float64, ReverseDiff.TrackedArray{Float64, Float64, 2, Matrix{Float64}, Matrix{Float64}}}})
Closest candidates are:
DiffEqOperators.BoundaryPaddedVector(::T, ::T, ::T2) where {T, T2<:AbstractVector{T}} at /Users/salazardetro1/.julia/packages/DiffEqOperators/FnmOt/src/boundary_padded_arrays.jl:14
```

I believe that my current equations might be solved using another method (maybe an ODEProblem with a mass matrix?), but my next step will be to include a nonlinear term multiplying the time derivatives as in

and I think that I can only solve those problems with a `DAEProblem`

. With that said, how can I compute the gradients? Thanks.