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,CairoMakie
    default_plotter!(CairoMakie)
    CairoMakie.activate!(type="svg")
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]
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])
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]
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)
    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=1:size(inival,2)
   	 		inival[1,i]=fpeak(coord[:,i])
   	 		inival[2,i]=0
            #
    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 = 12.953886341, nd = 1925, njac = 855, nf = 2780)