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 DataStructures
    using GridVisualize,CairoMakie
end

function barenblatt(x,t,m)
    tx=t^(-1.0/(m+1.0))
    xx=x*tx
    xx=xx*xx
    xx=1- xx*(m-1)/(2.0*m*(m+1));
    if xx<0.0
        xx=0.0
    end
    return tx*xx^(1.0/(m-1.0))
end

barenblatt (generic function with 1 method)

function create_porous_medium_problem(n,m)
    h=1.0/convert(Float64,n/2)
    X=collect(-1:h:1)
    grid=VoronoiFVM.Grid(X)

    function flux!(f,u,edge,data)
        f[1]=u[1,1]^m-u[1,2]^m
    end

    storage!(f,u,node,data)= f[1]=u[1]

    sys=VoronoiFVM.System(grid,flux=flux!,storage=storage!, species=1)
    sys,X
end

create_porous_medium_problem (generic function with 1 method)

begin
function run_vfvm(;n=20,m=2,t0=0.001, tend=0.01,tstep=1.0e-6)
    sys,X=create_porous_medium_problem(n,m)
    inival=unknowns(sys)
    inival[1,:].=map(x->barenblatt(x,t0,m),X)
    sol=VoronoiFVM.solve(sys;inival,times=(t0,tend),Δt=tstep,Δu_opt=0.01,Δt_min=tstep,store_all=true,log=true, reltol=1.0e-3)
    err=norm(sol[1,:,end]-map(x->barenblatt(x,tend,m),X))
    sol,sys,err
end
run_vfvm(m=2,n=10) # "Precompile"
end;

begin
    function run_diffeq(;n=20,m=2, t0=0.001,tend=0.01,solver=nothing)
        sys,X=create_porous_medium_problem(n,m)
        inival=unknowns(sys)
        inival[1,:].=map(x->barenblatt(x,t0,m),X)
        state=VoronoiFVM.SystemState(sys)
        problem = ODEProblem(state,inival,(t0,tend))
        odesol = solve(problem,solver)
        sol=reshape(odesol,sys; state)
        err=norm(sol[1,:,end]-map(x->barenblatt(x,tend,m),X))
        sol, sys,err
    end
    for method in diffeqmethods
        run_diffeq(m=2,n=10,solver=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

t1=@elapsed sol1,sys1,err1=run_vfvm(m=m,n=n);history_summary(sol1)

(seconds = 2.03, tasm = 1.67, tlinsolve = 0.338, steps = 832, iters = 1662, maxabsnorm = 2.65e-6, maxrelnorm = 0.000263, maxroundoff = 2.46e-16, iters_per_step = 2.0, facts_per_step = 0.0, liniters_per_step = 0.0)

m=2; n=50;

t2=@elapsed sol2,sys2,err2=run_diffeq(m=m,n=n,solver=diffeqmethods[method]());history_summary(sol2)

(nd = 165, njac = 82, nf = 248)

@test err2<err1

Test Passed