Got a working example! It seems to have something to do with this equation.

```
0.0 ~ 0.1 * u(t,x)*max(0,(u(t,x) - v(t,x)))^r - (u(t,x) - v(t,x)) * ((1.1+sin(t))^r)]
```

and in particular, this term `max(0,(u(t,x) - v(t,x)))^r`

```
using DifferentialEquations, ModelingToolkit, MethodOfLines, DomainSets, ForwardDiff, FiniteDiff, StatsBase
# Parameters, variables, and derivatives
@parameters t, x, p, q, r
@variables u(..), v(..)
#@variables u[1:n_comp](..)
Dt = Differential(t)
Dx = Differential(x)
Dxx = Differential(x)^2
params = Symbolics.scalarize(reduce(vcat,[p => 1.5, q => 1.2, r .=> 1.3]))
# 1D PDE and boundary conditions
eqs = [Dt(u(t, x)) ~ -q*Dx(u(t,x)) + p * Dxx(u(t, x)),
0.0 ~ 0.1 * u(t,x)*max(0,(u(t,x) - v(t,x)))^r - (u(t,x) - v(t,x)) * ((1.1+sin(t))^r)]
bcs = [u(0, x) ~ 0.0,
v(0,x) ~ 0.0,
u(t, 0) ~ 0.5*exp((-(t^2))/0.01),
v(t,0) ~ 0.5*exp((-(t^2))/0.01),
Dx(u(t,0)) ~ 0.0,
Dx(v(t,0)) ~ 0.0]
# Space and time domains
domains = [t ∈ Interval(0.0, 1.0),
x ∈ Interval(0.0, 1.0)]
# PDE system
@named pdesys = PDESystem(eqs, bcs, domains, [t,x], [u(t,x), v(t,x)], params)
dx = 0.1
# Method of lines discretization
order = 2
discretization = MOLFiniteDifference([x=>dx], t; approx_order = order)
# Convert the PDE problem into an ODE problem
prob = discretize(pdesys,discretization)
# Solve ODE problem
sol = solve(prob, KenCarp47(), saveat=0.02);
function lossfun(sol)
output = []
for i = 1:10
push!(output,mean(sol((i-1)/10:i/10, 1.0, dv = u(t,x))))
end
outputvec = stack(output)'
diffs = sum(abs.(outputvec.-0))
return diffs
end
function errortest(param_vec)
loss = []
newprob = remake(prob, p = param_vec)
try
sol = solve(newprob, saveat = 0.1, KenCarp47(), reltol=10^-8, abstol=10^-8)
if sol.retcode != :Success
println("anomalous solution generated")
loss = Inf
else
diffs = lossfun(sol)
loss = diffs
end
catch
println("solve failed")
push!(loss, Inf)
end
return loss
end
default_params = [1.5 1.2 1.3]
#This works
errortest(default_params)
#This does not
ForwardDiff.gradient(errortest,default_params)
#This works
FiniteDiff.finite_difference_gradient(errortest,default_params)
```