begin
    import Pkg as _Pkg
    haskey(ENV, "PLUTO_PROJECT") && _Pkg.activate(ENV["PLUTO_PROJECT"])
    using Test
    using Revise
    using Printf
    using VoronoiFVM
    using OrdinaryDiffEqBDF
    using OrdinaryDiffEqRosenbrock
    using OrdinaryDiffEqSDIRK
    using LinearAlgebra
    using PlutoUI
    using ExtendableGrids
    using DataStructures
    using GridVisualize
    if isdefined(Main, :PlutoRunner)
        using CairoMakie
        default_plotter!(CairoMakie)
        CairoMakie.activate!(type = "png")
    end
end

A Brusselator problem

begin
    const bruss_A = 2.25
    const bruss_B = 7.0
    const bruss_D_1 = 0.025
    const bruss_D_2 = 0.25
    const pert = 0.1
    const bruss_tend = 150
end;

function bruss_storage(f, u, node, data)
    f[1] = u[1]
    f[2] = u[2]
    return nothing
end;

function bruss_diffusion(f, u, edge, data)
    f[1] = bruss_D_1 * (u[1, 1] - u[1, 2])
    f[2] = bruss_D_2 * (u[2, 1] - u[2, 2])
    return nothing
end;

function bruss_reaction(f, u, node, data)
    f[1] = (bruss_B + 1.0) * u[1] - bruss_A - u[1]^2 * u[2]
    f[2] = u[1]^2 * u[2] - bruss_B * u[1]
    return nothing
end;

begin

    function ODESolver(system, inival, solver)
        state = VoronoiFVM.SystemState(system)
        problem = ODEProblem(state, inival, (0, bruss_tend))
        odesol = solve(
            problem,
            solver,
            dt = 1.0e-5, reltol = 1.0e-4
        )
        return reshape(odesol, system; state)
    end

    sys0 = VoronoiFVM.System(simplexgrid(0:0.1:1), species = [1, 2], flux = bruss_diffusion, storage = bruss_storage, reaction = bruss_reaction)
    problem0 = ODEProblem(sys0, unknowns(sys0), (0, 0.1))

    for method in diffeqmethods
        solve(problem0, method.second()) #precompile
    end
end

OrderedDict{String, UnionAll} with 4 entries:
  "Rosenbrock23 (Rosenbrock)"    => Rosenbrock23
  "QNDF2 (Like matlab's ode15s)" => QNDF2
  "FBDF"                         => FBDF
  "Implicit Euler"               => ImplicitEuler

if bruss_dim == 1
    bruss_X = -1:0.01:1
    bruss_grid = simplexgrid(bruss_X)
else
    bruss_X = -1:0.1:1
    bruss_grid = simplexgrid(bruss_X, bruss_X)
end;

bruss_system = VoronoiFVM.System(
    bruss_grid, species = [1, 2],
    flux = bruss_diffusion, storage = bruss_storage, reaction = bruss_reaction
);

begin
    inival = unknowns(bruss_system, inival = 0)
    coord = bruss_grid[Coordinates]
    fpeak(x) = exp(-norm(10 * x)^2)
    for i in 1:size(inival, 2)
        @views inival[1, i] = fpeak(coord[:, i])
        @views inival[2, i] = fpeak(coord[:, i])
    end
end

t_run = @elapsed bruss_tsol = ODESolver(bruss_system, inival, diffeqmethods[bruss_method]());

(t_run = t_run, VoronoiFVM.history_details(bruss_tsol)...)

(t_run = 5.14897588, nd = 1775, njac = 787, nf = 2562)