Dear all, I have run into the following issue, which is described in details also here: Frequently Asked Questions · DifferentialEquations.jl (sciml.ai).

I have tried all suggestions on the link and made a great amount of small modifications to test alternatives. Also, the same exact implementation presents no instability in MATHEMATICA, so I strongly believe the model is correct.

```
Warning: dt <= dtmin. Aborting. There is either an error in your model specification or the true solution is unstable.
└ @ SciMLBase C:\Users\18143\.julia\packages\SciMLBase\L7Nun\src\integrator_interface.jl:345
```

Here is the main part of my code relevant to this:

```
##
##Define the parameters of the dynamics
dt_step=0.10
time=collect(dt_step:dt_step:4.0)
initialcond=collect(0.1:0.2:0.9)
N=4; M=4; T=4; L=1; Nsamples=(N+1)*(M+1)
##Problem-related parameters
β=0.5; ξ=0.5; α=-2; C=2;
### --- Functions --- ###
function uₜₑᵣₘₛ(i,j,t,x₀)
uₜₑᵣₘₛ = aₘₙ[i,j]*cos((i-1)*π*x₀/L)*cos((j-1)*π*t/T)
end
function uₐₚₚ(t,x₀)
uₐₚₚ = sum(uₜₑᵣₘₛ(i,j,t,x₀) for i=1:(M+1), j=1:(N+1))
end
### --- More Functions --- ###
using DifferentialEquations
function ϕᵣₑₐₗ(x₀)
f(ϕ,p,t) = β*ϕ*(1-ϕ)*(uₐₚₚ(t,x₀)-ξ)
ϕ₀=x₀
tspan=(0.0,4.0)
prob=ODEProblem(f,ϕ₀,tspan)
sol=solve(prob, Tsit5(), reltol=1e-8, abstol=1e-8) #Tested many of these already
return sol
end
#Numerical Integration of Functional - Trapezoidal Rule
function quad_trap(f,a,b,NN,x₀)
h = (b-a)/NN
int = h * ( f(a,x₀) + f(b,x₀) ) / 2
for k=1:NN-1
xk = (b-a) * k/NN + a
int = int + h*f(xk,x₀)
end
return int
end
function RealFunctional(t,x₀)
α*ϕᵣₑₐₗ(x₀)(t)*uₐₚₚ(t,x₀) + C*uₐₚₚ(t,x₀)^2
end
function RealFunctionalCost(x₀)
Int = quad_trap(RealFunctional,0,T,size(time)[1],x₀)
return Int
end
function RealJₐ()
return sum(RealFunctionalCost(x₀) for x₀=initialcond)
end
#Gradient Descent
aₘₙ=rand(Uniform(-1,1),(M+1),(N+1)); ϵ=10e-2;
ApproxRealGrad=Matrix{Float64}(undef,(M+1),(N+1))
#Gradient Computation -- I know this is lousy but it needs to be like this
FunctionalI=RealJₐ()
println(FunctionalI)
a_δ=aₘₙ
for i=1:(M+1), j=1:(N+1)
aₘₙ=a_δ
aₘₙ[i,j] = a_δ[i,j] + ϵ
FunctionalF=RealJₐ()
ApproxRealGrad[i,j] = (FunctionalF-FunctionalI)/ϵ
end
aₘₙ=a_δ
τ=0.1
#Descent
aₘₙ -= τ*ApproxRealGrad
```

There is a loop in the above omitted here, because I am running iteration by iteration “by hand” in Jupyter.

Also, I know there would be many alternatives to the computation of this derivative, but for our application this is a necessary step. The same step in the descent step worked fine in MATHEMATICA as well.

Please help me, mainly with the ODE `ϕᵣₑₐₗ(x₀)`

, which is the source of this.

I spent a lot of time already trying to get out of this warning. Thank you a lot for your time and help.

Gabriel