Hello everyone, thank you for all your advice. Now I have a solution for this type of problem. I will just give the code, and refer to this paper for basic knowledge about the Fillippov systems.

```
using StaticArrays, OrdinaryDiffEq, Plots, LinearAlgebra
struct ODESol
solution # ODE solution solved by OrdinaryDiffEq.jl
state::Symbol # only two symbols are used: :regularmode and :slidemode
end
function fillippov_solve(f1, f2, h, ∇h, v0, p0, t0, tend)
# it is assumed that f1 is at h(x)<0 and f2 is at h(x)>0
# f1, f2, h, ∇h are all type of f(u, p, t), the output is a vector
# initial condition: v0
# time interval: t0,tend
# parameter p0
#
# first define the vector field on the switch surface
function sv(x, p, t)
α = dot(∇h(x, p, t), f1(x, p, t)) / dot(∇h(x, p, t), f1(x, p, t) - f2(x, p, t))
(1 - α) * f1(x, p, t) + α * f2(x, p, t)
end
# the hyper surface
function newh(x, t, integrator)
h(x, integrator.p, t)
end
# stop the integrator
function stop(integrator)
if (integrator.t-integrator.sol.t[1])< 1e-2 # to avoid the integrator stopping at the beginning when it exits the sliding surface
nothing
else
terminate!(integrator)
end
end
# exit the sliding mode when slide_exit=0
function slide_exit(x, t, integrator)
dot(∇h(x, integrator.p, t), f2(x, integrator.p, t)) * dot(∇h(x, integrator.p, t), f1(x, integrator.p, t))
end
# main data consist of different solutions at h>0, h<0 and h=0
sols = Vector{ODESol}(undef, 1)
if abs(h(v0, p0, t0)) < 1e-10 # if this happen, we assume that the initial condition is on h=0
prob = ODEProblem(sv, v0, (t0, tend), p0)
cb = ContinuousCallback(slide_exit, stop)
firstsol = ODESol(solve(prob, Vern6(), callback=cb), :slidemode)
sols[1] = firstsol
elseif h(v0, p0, t0) > 0
prob = ODEProblem(f2, v0, (t0, tend), p0)
cb = ContinuousCallback(newh, stop)
firstsol = ODESol(solve(prob, Vern6(), callback=cb), :regularmode)
sols[1] = firstsol
else
prob = ODEProblem(f1, v0, (t0, tend), p0)
cb = ContinuousCallback(newh, stop)
firstsol = ODESol(solve(prob, Vern6(), callback=cb), :regularmode)
sols[1] = firstsol
end
while sols[end].solution.t[end] != tend # if the last of sols's time ≠ tend
v1 = sols[end].solution.u[end]
t1 = sols[end].solution.t[end]
if sols[end].state == :slidemode # if the last solution is in slidemode, we have to determine which vector field it enter
if abs(dot(∇h(v1, p0, t1), f2(v1, p0, t1))) < 1e-8 && dot(∇h(v1, p0, t1), f1(v1, p0, t1)) > 1e-5
prob = ODEProblem(f2, v1, (t1, tend), p0)
cb = ContinuousCallback(newh, stop, repeat_nudge=0.2)
firstsol = ODESol(solve(prob, Vern6(), callback=cb), :regularmode)
append!(sols, [firstsol])
elseif abs(dot(∇h(v1, p0, t1), f1(v1, p0, t1))) < 1e-8 && dot(∇h(v1, p0, t1), f2(v1, p0, t1)) < -1e-5
prob = ODEProblem(f1, v1, (t1, tend), p0)
cb = ContinuousCallback(newh, stop, repeat_nudge=0.2)
firstsol = ODESol(solve(prob, Vern6(), callback=cb), :regularmode)
append!(sols, [firstsol])
else
break
print("The solution is not defined on state $v1 and time $t1")
end
else
if dot(∇h(v1, p0, t1), f1(v1, p0, t1)) * dot(∇h(v1, p0, t1), f2(v1, p0, t1)) > 0 # crossing h=0
if dot(∇h(v1, p0, t1), f1(v1, p0, t1)) > 0
prob = ODEProblem(f2, v1, (t1, tend), p0)
cb = ContinuousCallback(newh, stop, repeat_nudge=0.2)
firstsol = ODESol(solve(prob, Vern6(), callback=cb), :regularmode)
append!(sols, [firstsol])
else
prob = ODEProblem(f1, v1, (t1, tend), p0)
cb = ContinuousCallback(newh, stop, repeat_nudge=0.2)
firstsol = ODESol(solve(prob, Vern6(), callback=cb), :regularmode)
append!(sols, [firstsol])
end
else # entering h=0
prob = ODEProblem(sv, v1, (t1, tend), p0)
cb = ContinuousCallback(slide_exit, stop)
firstsol = ODESol(solve(prob, Vern6(), callback=cb), :slidemode)
append!(sols, [firstsol])
end
end
end
sols
end
```

The main function is `fillippov_solve`

, and its output is a vector consist of `ODESol`

.

Now we can test this function using the system I gave before. First, define two vector fields:

```
function F1(x)
A, γ, λ, η, α = 10.0, 3.0, 0.1, 0.01, 0.3
A * (α / (1 + γ * abs(x)) + λ + η * x^2)
end
function F2(x)
A, γ, λ, η, α = 10.0, 3.0, 0.1, 0.01, 0.3
-A * (α / (1 + γ * abs(x)) + λ + η * x^2)
end
function vec1(x, p, t)
SVector(x[2], -p[1] * x[2] - p[2] * x[1] + F1(x[2] - p[3]) + p[4] * cos(p[5] * t)) # p=[c,k,vs,ϵ,ω]
end
function vec2(x, p, t)
SVector(x[2], -p[1] * x[2] - p[2] * x[1] + F2(x[2] - p[3]) + p[4] * cos(p[5] * t)) # p=[c,k,vs,ϵ,ω]
end
```

And then define the hypersurface and its grad:

```
function H(x, p, t)
x[2] - p[3]
end
function ∇H(x, p, t)
SVector(0.0, 1.0)
end
```

Note that we have to make sure that when H>0 the vector field is `vec2`

and when H<0 the vector field is `vec1`

.

Now we can use `fillippov_solve`

to solve this ODE. Just run

```
sol = fillippov_solve(vec1, vec2, H, ∇H, SVector(1.0, 3.0), SVector(0.2, 1, 1, 0.2, 1.067), 0.0, 20.0)
```

You will get a vector.

```
sol[2].state
```

gives

```
:slidemode
```

```
sol[2].solution
```

gives a ODE solution

```
retcode: Terminated
Interpolation: specialized
6th order lazy interpolation
t: 3-element Vector{Float64}:
0.5146478247426896
2.117126414261962
2.117126414261962
u: 3-element Vector{StaticArraysCore.SVector{2, Float64}}:
[2.07049529563982, 1.0000000000000895]
[3.6729738851570852, 1.000000000000007]
[3.6729738851570852, 1.000000000000007]
```

Since we define a new vector field `sv`

on the switch surface when the vector fields are squeezing towards the switch surface, the solution is now stable.

We can define a simple function to plot all the solutions:

```
function myplot(sols)
plotlyjs()
figure = plot()
for i in sols
plot!(i.solution, idxs=(1,2))
end
figure
end
```

Now we run

```
myplot(sol)
```

it gives this picture:

Since I also consider the case of crossing, the codes should work on the usual piecewise smooth case.