using IVPTestSuite  # imports ODE, DASSL and Sundials
const S = Solvers # for convenience

# choose test-case (all of them are in IVPTestSuite.tc_all)
tc_name = [:vdpol
           :rober
           :threebody
           :bruss1d
           :chemakzo
           :hires][3]
tc = IVPTestSuite.tc_all[tc_name]           

# Pick a solver.  S.allsolvers is the list of all solvers.  Generally
# each family of solvers also has a list.
solver = S.ode45_dp

###
# make a TestRun which combines a TestCase with a Solver + some extras:
abstol = 1e-7
reltol = abstol
dt0 = NaN # size of first step
tr = TestRunAdapt(tc, solver, Dict{Symbol,Any}(), abstol, reltol, dt0)
# For fixed step solvers use TestRunFixedStep instead.

# run it and display results:
re = IVPTestSuite.run_ode_test_throwerror(tr) # this returns a TestResults instance
re = IVPTestSuite.run_ode_test_throwerror(tr) # Run it twice to get accurate timings
show(re)

using IVPTestSuite  # imports ODE, DASSL and Sundials
const S = Solvers # for convenience

# choose test-case (all of them are in IVPTestSuite.tc_all)
tc_name = [:vdpol
           :rober
           :threebody
           :bruss1d
           :chemakzo
           :hires][3]
tc = IVPTestSuite.tc_all[tc_name]           

# Pick a solver.  S.allsolvers is the list of all solvers.  Generally
# each family of solvers also has a list.
solver = S.ode45_dp

###
# make a TestRun which combines a TestCase with a Solver + some extras:
abstols = 10.0.^(-5:-1:-10)
reltols = abstols
dt0 = NaN # size of first step
suite = TestSuite(tc, solver, abstols, reltols, [NaN])

# running it will display results:
res = run_ode_testsuite(suite)

# create an instance of TestCaseExplicit:
mytest = let
    tcname = ...::Symbol # name should be same as variable name (except for upper/lower case)
    
    T = ... # the datatype used, probably Float64
    Tarr = [Matrix, SparseMatrixCSC][1] # the container datatype used
                                        # for the Jacobian and mass
                                        # matrix, probably Matrix
                                        # (default) or
                                        # SparseMatrixCSC.
    dof = ...::Int # degrees of freedom
    dae = ...::Int  # index of DAE, ==0 for ODE
    # stiffness of system, one of the three constants
    stiffness = [nonstiff, mildlystiff, stiff][...]

    ## the problem function
    function fn!(t::T,y::Vector{T},dydt::Vector{T})
        # The ode function dydt = f(t,y) modifying dydt in-place.
        # Note that all indices of dydt need to be written!
        ...
        return nothing
    end
    # initializes storage for y:
    fn!( ; T_::Type=T, dof_=dof) = zeros(T_,dof_)
 
    # NOTE, using keyword arguments more intuitively like:
    # fn!( ; T::Type=T, dof=dof) = zeros(T,dof)
    # does not work: https://github.com/JuliaLang/julia/issues/9948

    ## Jacobian of fn
    function jac!(t::T,y::Vector{T},dfdy::Tarr)
        # the Jacobian of fn, inplace
        #
        # Set jac=nothing if it is not known.
        #
        # Note that jac! should write all indices of given in
        # jpattern.
        ...
        return nothing
    end
    # Returns a matrix which can hold the Jacobian of selected
    # type, initialized to zero.  Can also be used to make a
    # Matlab-style JPattern matrix with jac!(Bool)
    jac!( ; T_::Type=T, dof_=dof) = zeros(T_,dof_,dof_)  # if the problem is large better return an appropriate sparse matrix

    ## Mass matrix
    function mass!(t::T,y::Vector{T},m::Tarr)
        # mass matrix:  M(t,y) * dydt = f(t,y)
        #
        # Set to nothing if it is I, i.e. mass=nothing
        ...
        return nothing
    end
    # Returns a matrix which can hold the mass matrix.
    # Can also be used to make a Matlab-style MvPattern matrix with mass!(Bool)
    mass!( ; T_::Type=T, dof_=dof) = zeros(T_,dof_,dof_)  # if the problem is large better return an appropriate sparse matrix
    
    ic = T[...] # vector of initial conditions
    tspan = T[start, stop] # integration interval
    refsol = T[...] # reference solution at tspan[2]

    refsolinds = trues(dof)   # if refsol does not contain all indices
                              # then specify which indices of the
                              # solution are to be compared to refsol
    
    scd_absinds = Int[]       # set where refsol is very small or
                              # zero. To avoid it dominating the
                              # relative error
    
    tc = TestCaseExplicit{tcname, T, Tarr}(
                             stiffness,
                             dae,
                             dof,
                             fn!,
                             jac!,
                             mass!,
                             ic,
                             tspan,
                             refsol,
                             refsolinds,
                             scd_absinds)

    # put into the right buckets
    tc_all[tcname] = tc
    tc_stiff[tcname] = tc
    # Note that the return-result of the last statement of the let
    # block needs to be tc!
    tc
end

# Fill in the ...
function wrapped_solver(tr::TestRun)
    # Wraps the specific MyPkg.mysolver such that it works within
    # IVPTestSuite setup.
    tc = tr.tc
    so = tr.solver
    ###
    # 0) Wrap tc.fn!, tc.jac!, tc.mass! if necessary
    
    ###
    # 1) Make call signature

    args = ...
    kwargs = ...
    
    ###
    # 2) Call solver, if it does not succeed throw an error (if that
    # is not done anyway)
    out = so.solverfn(args...; kwargs...)  
    # (probably no need to modify this section)
    
    ###
    # 3) Transform output to conform to standard:
    # tend -- end time reached
    # yend -- solution at tend
    # stats -- statistics, if available: (steps_total,steps_accepted, fn_evals, jac_evals, linear_solves)
    #                         otherwise  (-1, -1, -1, -1, -1)
    ...
    return tend, yend, stats
end
    
solverfn = ...::Function # the actual solver function, for example ODE.ode23s.
# choose from 
typ =  [:ex, :im, :imex][2]                    # Typ: :ex (explicit method), :im (implicit method), :imex (IMEX method)
stiffness = [nonstiff, mildlystiff, stiff][1]  # stiff, mildly stiff or non-stiff solver
adaptive = [true, false][1]                    # is the solver adaptive
daeindex = 0                                   # maximum index of DAE it can solve, ==0 for ODE only
eq_type = [ explicit_eq, explicit_mass_eq, implicit_eq][1] # [dy/dt=F(t,y); M(t,y)*dy/dt=F(t,y); F(t,y,dy/dt)=0]

solverfn = MyPkg.myode
solverpkg = MyPkg

myodesolver = Solver{typ}(solverfn, solverpkg, wrapped_solver, stiffness, adaptive, daeindex, explicit_eq)

push!(allsolvers, myodesolver)