Fitzhugh-Nagumo Work-precision Diagrams

using Pkg, OrdinaryDiffEq, ParameterizedFunctions, DiffEqDevTools, Plots

f = @ode_def FitzhughNagumo begin
  dv = v - v^3/3 -w + l
  dw = τinv*(v +  a - b*w)
end a b τinv l

p = [0.7,0.8,1/12.5,0.5]
prob = ODEProblem(f,[1.0;1.0],(0.0,10.0),p)

abstols = 1.0 ./ 10.0 .^ (6:13)
reltols = 1.0 ./ 10.0 .^ (3:10);

sol = solve(prob,Vern7(),abstol=1/10^14,reltol=1/10^14)
test_sol = TestSolution(sol)
using Plots; gr()

plot(sol,dpi=200,linewidth=1,legend=:bottomleft,legendfontsize=6)

setups = [
          Dict(:alg=>DP5()),
          Dict(:alg=>BS5()),
          Dict(:alg=>Tsit5()),
          Dict(:alg=>Vern6()),
          #Dict(:alg=>RKC()),
          Dict(:alg=>ROCK2()),
          Dict(:alg=>ROCK4()),
          Dict(:alg=>SERK2v2(controller=:PI)),
          Dict(:alg=>SERK2v2(controller=:Predictive)),
          Dict(:alg=>ESERK5())
          ]
#Names = ["DP5" "BS5" "Tsit5" "Vern6" "RKC" "ROCK2" "ROCK4" "SERK2 PI" "SERK2 Predictive" "ESERK5"]
Names = ["DP5" "BS5" "Tsit5" "Vern6" "ROCK2" "ROCK4" "SERK2 PI" "SERK2 Predictive" "ESERK5"]

wp = WorkPrecisionSet(prob,abstols,reltols,setups;names=Names,appxsol=test_sol,save_everystep=false,numruns=100,maxiters=Int(1e6))
plot(wp,dpi=200,linewidth=1,legend=:topright,legendfontsize=6)

wp = WorkPrecisionSet(prob,abstols,reltols,setups;names=Names,appxsol=test_sol,save_everystep=false,numruns=100,maxiters=Int(1e6),error_estimate=:L2,dense_errors=true)
plot(wp,dpi=200,linewidth=1,legend=:topright,legendfontsize=6)

setups = [
          Dict(:alg=>DP5()),
          #Dict(:alg=>BS5()),
          Dict(:alg=>Tsit5()),
          Dict(:alg=>Vern6()),
          #Dict(:alg=>RKC()),
          Dict(:alg=>ROCK2()),
          #Dict(:alg=>ROCK4()),
          Dict(:alg=>SERK2v2(controller=:PI)),
          Dict(:alg=>SERK2v2(controller=:Predictive)),
          Dict(:alg=>ESERK5())
          ]
Names = ["DP5" "Tsit5" "Vern6" "ROCK2" "SERK2 PI" "SERK2 Predictive" "ESERK5"]
wp = WorkPrecisionSet(prob,abstols,reltols,setups;names=Names,appxsol=test_sol,save_everystep=false,numruns=100,maxiters=Int(1e6),error_estimate=:L2,dense_errors=true)
plot(wp,dpi=200,linewidth=1,legend=:topleft,legendfontsize=6)