using ReachabilityAnalysis, Plots

@taylorize function seir2!(du,u,p,t)
  S, E, I, R, α, β, γ = u
  βIS = β * (I * S)
  αE = α*E
  γI = γ*I
  du[1] = -βIS      # dS
  du[2] = βIS - αE  # dE
  du[3] = -γI + αE  # dI
  du[4] = γI        # dR
  local zerou = zero(u[1])
  du[5] = zerou
  du[6] = zerou
  du[7] = zerou
end

E₀ = 1e-4
u₀ = [1-E₀, E₀, 0, 0]
param_error = 0.01
α = 0.2 ± param_error
γ = 0.5 ± param_error


Threads.@threads for i = 1:10
    β = i
    p = [α, β, γ]
    X0 = IntervalBox(vcat(u₀, p));
    prob = @ivp(x' = seir2!(x), dim=7, x(0) ∈ X0);
    sol = solve(prob, tspan=(0.0, 200.0), alg=TMJets21a(orderT=7, orderQ=1));
end

.............$ julia -t 4 test.jl
ERROR: LoadError: TaskFailedException
Stacktrace:
 [1] wait
   @ ./task.jl:334 [inlined]
 [2] threading_run(func::Function)
   @ Base.Threads ./threadingconstructs.jl:38
 [3] top-level scope
   @ ./threadingconstructs.jl:97

    nested task error: AssertionError: order ≤ get_order() && lencoef ≤ num_coeffs
    Stacktrace:
      [1] resize_coeffsHP!(coeffs::Vector{Float64}, order::Int64)
        @ TaylorSeries ~/.julia/packages/TaylorSeries/vTIIX/src/auxiliary.jl:39
      [2] HomogeneousPolynomial
        @ ~/.julia/packages/TaylorSeries/vTIIX/src/constructors.jl:100 [inlined]
      [3] HomogeneousPolynomial
        @ ~/.julia/packages/TaylorSeries/vTIIX/src/constructors.jl:106 [inlined]
      [4] HomogeneousPolynomial(#unused#::Type{Float64}, nv::Int64)
        @ TaylorSeries ~/.julia/packages/TaylorSeries/vTIIX/src/constructors.jl:132
      [5] #TaylorN#27
        @ ~/.julia/packages/TaylorSeries/vTIIX/src/constructors.jl:202 [inlined]
      [6] TaylorN
        @ ~/.julia/packages/TaylorSeries/vTIIX/src/constructors.jl:202 [inlined]
      [7] set_variables(::Type{Float64}, names::Vector{String}; order::Int64)
        @ TaylorSeries ~/.julia/packages/TaylorSeries/vTIIX/src/parameters.jl:146
      [8] set_variables(::Type{Float64}, names::String; order::Int64, numvars::Int64)
        @ TaylorSeries ~/.julia/packages/TaylorSeries/vTIIX/src/parameters.jl:166
      [9] #set_variables#15
        @ ~/.julia/packages/TaylorSeries/vTIIX/src/parameters.jl:171 [inlined]
     [10] post(alg::TMJets21a{Float64, ZonotopeEnclosure}, ivp::InitialValueProblem{BlackBoxContinuousSystem{typeof(seir2!)}, IntervalBox{7, Float64}}, timespan::IntervalArithmetic.Interval{Float64}; Δt0::IntervalArithmetic.Interval{Float64}, kwargs::Base.Pairs{Symbol, Any, Tuple{Symbol, Symbol}, NamedTuple{(:tspan, :alg), Tuple{Tuple{Float64, Float64}, TMJets21a{Float64, ZonotopeEnclosure}}}})
        @ ReachabilityAnalysis ~/.julia/packages/ReachabilityAnalysis/HmeX4/src/Algorithms/TMJets/TMJets21a/post.jl:29
     [11] solve(::InitialValueProblem{BlackBoxContinuousSystem{typeof(seir2!)}, IntervalBox{7, Float64}}; kwargs::Base.Pairs{Symbol, Any, Tuple{Symbol, Symbol}, NamedTuple{(:tspan, :alg), Tuple{Tuple{Float64, Float64}, TMJets21a{Float64, ZonotopeEnclosure}}}})
        @ ReachabilityAnalysis ~/.julia/packages/ReachabilityAnalysis/HmeX4/src/Continuous/solve.jl:64
     [12] macro expansion
        @ ~/YXY/8.0/case4/test.jl:30 [inlined]
     [13] (::var"#51#threadsfor_fun#1"{UnitRange{Int64}})(onethread::Bool)
        @ Main ./threadingconstructs.jl:85
     [14] (::var"#51#threadsfor_fun#1"{UnitRange{Int64}})()
        @ Main ./threadingconstructs.jl:52
in expression starting at /home/dso/YXY/8.0/case4/test.jl:25