diff --git a/News.md b/News.md
index cc1a5882..7082d627 100644
--- a/News.md
+++ b/News.md
@@ -3,6 +3,10 @@ BifurcationKit.jl, Changelog
 
 All notable changes to this project will be documented in this file (hopefully). No performance improvements will be notified but mainly the addition of new methods, the modifications of internal structs, etc.
 
+## [0.4.2]
+- change bordered linear solvers' interface
+- `record_from_solution` has been changed to the following definition
+
 ## [0.4.0]
 - Setfield.jl is not anymore the main component of BifurcationKit to support parameter axis. It has been changed in favour of Accessors.jl
 
diff --git a/examples/COModel.jl b/examples/COModel.jl
index 1991d129..fd664e55 100644
--- a/examples/COModel.jl
+++ b/examples/COModel.jl
@@ -19,7 +19,7 @@ par_com = (q1 = 2.5, q2 = 2.0, q3 = 10., q4 = 0.0675, q5 = 1., q6 = 0.1, k = 0.4
 
 z0 = [0.001137, 0.891483, 0.062345]
 
-prob = BifurcationProblem(COm!, z0, par_com, (@optic _.q2); record_from_solution = (x, p) -> (x = x[1], y = x[2], s = x[3]))
+prob = BifurcationProblem(COm!, z0, par_com, (@optic _.q2); record_from_solution = (x, p; k...) -> (x = x[1], y = x[2], s = x[3]))
 
 opts_br = ContinuationPar(dsmax = 0.015, dsmin=1e-4, ds=1e-4, p_min = 0.5, p_max = 2.0, n_inversion = 6, detect_bifurcation = 3, nev = 3)
 br = @time continuation(prob, PALC(), opts_br;
@@ -27,7 +27,7 @@ br = @time continuation(prob, PALC(), opts_br;
     normC = norminf,
     bothside = true)
 
-BK.plot(br)#markersize=4, legend=:topright, ylims=(0,0.16))
+plot(br)#markersize=4, legend=:topright, ylims=(0,0.16))
 ####################################################################################################
 # periodic orbits
 function plotSolution(x, p; k...)
@@ -49,7 +49,7 @@ function plotSolution(ax, x, p; ax1 = nothing, k...)
     scatter!(ax, xtt.t, xtt[1,:]; markersize = 1.5, k...)
 end
 
-args_po = (    record_from_solution = (x, p) -> begin
+args_po = (    record_from_solution = (x, p; k...) -> begin
         xtt = BK.get_periodic_orbit(p.prob, x, @set par_com.q2 = p.p)
         return (max = maximum(xtt[1,:]),
                 min = minimum(xtt[1,:]),
@@ -95,6 +95,7 @@ brpo = @time continuation(br, 5, opts_po_cont,
     # alg = PALC(),
     # alg = MoorePenrose(tangent=PALC(tangent = Bordered()), method = BK.direct),
     δp = 0.0005,
+    linear_algo = COPBLS(),
     callback_newton = callbackCO,
     args_po...
     )
@@ -114,8 +115,8 @@ sn_codim2 = continuation(br, 3, (@optic _.k), ContinuationPar(opts_br, p_max = 3
     bdlinsolver = MatrixBLS()
     )
 
-BK.plot(sn_codim2)#, real.(sn_codim2.BT), ylims = (-1,1), xlims=(0,2))
-BK.plot(sn_codim2, vars=(:q2, :x), branchlabel = "Fold", plotstability = false);plot!(br,xlims=(0.8,1.8))
+plot(sn_codim2)#, real.(sn_codim2.BT), ylims = (-1,1), xlims=(0,2))
+plot(sn_codim2, vars=(:q2, :x), branchlabel = "Fold", plotstability = false);plot!(br,xlims=(0.8,1.8))
 
 hp_codim2 = continuation((@set br.alg.tangent = Bordered()), 2, (@optic _.k), ContinuationPar(opts_br, p_min = 0., p_max = 2.8, detect_bifurcation = 0, ds = -0.0001, dsmax = 0.08, dsmin = 1e-4, n_inversion = 6, detect_event = 2, detect_loop = true, max_steps = 50, detect_fold=false) ; plot = true,
     verbosity = 0,
@@ -127,6 +128,6 @@ hp_codim2 = continuation((@set br.alg.tangent = Bordered()), 2, (@optic _.k), Co
     bothside = true,
     bdlinsolver = MatrixBLS())
 
-BK.plot(sn_codim2, vars=(:q2, :x), branchlabel = "Fold", plotcirclesbif = true)
+plot(sn_codim2, vars=(:q2, :x), branchlabel = "Fold", plotcirclesbif = true)
 plot!(hp_codim2, vars=(:q2, :x), branchlabel = "Hopf",plotcirclesbif = true)
 plot!(br,xlims=(0.6,1.5))
\ No newline at end of file
diff --git a/examples/SH2d-fronts-cuda.jl b/examples/SH2d-fronts-cuda.jl
index ac9f98ef..cf4a835a 100644
--- a/examples/SH2d-fronts-cuda.jl
+++ b/examples/SH2d-fronts-cuda.jl
@@ -115,7 +115,7 @@ par = (l = -0.1, ν = 1.3, L = L)
 prob = BK.BifurcationProblem(F_shfft, AF(sol0), par, (@optic _.l) ;
     J =  J_shfft,
     plot_solution = (x, p;kwargs...) -> plotsol!(x; color=:viridis, kwargs...),
-    record_from_solution = (x, p) -> norm(x))
+    record_from_solution = (x, p; k...) -> norm(x))
 
 opt_new = NewtonPar(verbose = true, tol = 1e-6, linsolver = L, eigsolver = Leig)
 sol_hexa = @time newton(prob, opt_new, normN = norminf);
diff --git a/examples/SH2d-fronts.jl b/examples/SH2d-fronts.jl
index 8d5755d8..237bf826 100644
--- a/examples/SH2d-fronts.jl
+++ b/examples/SH2d-fronts.jl
@@ -57,7 +57,7 @@ par = (l = -0.1, ν = 1.3, L1 = L1);
 
 optnew = NewtonPar(verbose = true, tol = 1e-8, max_iterations = 20)
 
-prob = BifurcationProblem(F_sh, vec(sol0), par, (@optic _.l); J = dF_sh, plot_solution = (x, p; kwargs...) -> (plotsol!((x); label="", kwargs...)),record_from_solution = (x, p) -> (n2 = norm(x), n8 = norm(x, 8)), d2F=d2F_sh, d3F=d3F_sh)
+prob = BifurcationProblem(F_sh, vec(sol0), par, (@optic _.l); J = dF_sh, plot_solution = (x, p; kwargs...) -> (plotsol!((x); label="", kwargs...)),record_from_solution = (x, p; k...) -> (n2 = norm(x), n8 = norm(x, 8)), d2F=d2F_sh, d3F=d3F_sh)
 # optnew = NewtonPar(verbose = true, tol = 1e-8, max_iterations = 20, eigsolver = EigArpack(0.5, :LM))
 sol_hexa = @time newton(prob, optnew)
 println("--> norm(sol) = ", norm(sol_hexa.u, Inf64))
@@ -97,7 +97,7 @@ br = @time continuation(
 ###################################################################################################
 # codim2 Fold continuation
 optfold = @set optcont.detect_bifurcation = 0
-@set! optfold.newton_options.verbose = true
+@reset optfold.newton_options.verbose = true
 optfold = setproperties(optfold; p_min = -2., p_max= 2., dsmax = 0.1)
 
 # dispatch plot to fold solution
@@ -127,7 +127,7 @@ function dF_sh2(du, u, p)
 end
 
 prob2 = @set prob.VF.J = (u, p) -> (du -> dF_sh2(du, u, p))
-@set! prob2.u0 = vec(sol0)
+@reset prob2.u0 = vec(sol0)
 
 sol_hexa = @time newton(prob2, @set optnew.linsolver = ls)
 println("--> norm(sol) = ", norm(sol_hexa.u, Inf64))
@@ -168,7 +168,11 @@ plot!(br)
 deflationOp = DeflationOperator(2, 1.0, [sol_hexa.u])
 optcontdf = @set optcont.newton_options.verbose = false
 
-algdc = BK.DefCont(deflation_operator = deflationOp, perturb_solution = (x,p,id) -> (x  .+ 0.1 .* rand(length(x))), max_iter_defop = 50, max_branches = 40, seek_every_step = 5,)
+algdc = BK.DefCont(deflation_operator = deflationOp, 
+            perturb_solution = (x,p,id) -> (x  .+ 0.1 .* rand(length(x))), 
+            max_iter_defop = 50, 
+            max_branches = 40, 
+            seek_every_step = 5,)
 
 brdf = continuation(prob, algdc, setproperties(optcontdf; detect_bifurcation = 0, plot_every_step = 1);
     plot = true, verbosity = 2,
diff --git a/examples/SH3d.jl b/examples/SH3d.jl
index f7663de4..af8801cf 100644
--- a/examples/SH3d.jl
+++ b/examples/SH3d.jl
@@ -117,12 +117,12 @@ prob = BK.BifurcationProblem(F_sh, AF(vec(sol0)), par, (@optic _.l),
     J = (x, p) -> (dx -> dF_sh(x, p, dx)),
     # J = (x, p) -> J_sh(x, p),
     plot_solution = (ax, x, p; ax1=nothing) -> contour3dMakie!(ax, x),
-    record_from_solution = (x, p) -> (n2 = norm(x), n8 = norm(x, 8)),
+    record_from_solution = (x, p; k...) -> (n2 = norm(x), n8 = norm(x, 8)),
     issymmetric = true)
 
 optnew = NewtonPar(verbose = true, tol = 1e-8, max_iterations = 20, linsolver = @set ls.verbose = 0)
-@set! optnew.eigsolver = eigSH3d
-# @set! optnew.linsolver = DefaultLS()
+@reset optnew.eigsolver = eigSH3d
+# @reset optnew.linsolver = DefaultLS()
 sol_hexa = @time newton(prob, optnew)
 println("--> norm(sol) = ", norm(sol_hexa.u, Inf64))
 
@@ -173,11 +173,11 @@ br = @time continuation(
     # end
     )
 
-plot(br)
+BK.plot(br)
 ####################################################################################################
 get_normal_form(br, 3; nev = 25)
 
-br1 = @time continuation(br, 3, setproperties(optcont; save_sol_every_step = 10, detect_bifurcation = 3, p_max = 0.1, plot_every_step = 5, dsmax = 0.01);
+br1 = @time continuation(br, 3, setproperties(optcont; save_sol_every_step = 10, detect_bifurcation = 0, p_max = 0.1, plot_every_step = 5, dsmax = 0.01);
     plot = true, verbosity = 3,
     δp = 0.005,
     verbosedeflation = true,
diff --git a/examples/SHpde_snaking.jl b/examples/SHpde_snaking.jl
index 61880dcd..a91053c6 100644
--- a/examples/SHpde_snaking.jl
+++ b/examples/SHpde_snaking.jl
@@ -30,7 +30,7 @@ parSH = (λ = -0.1, ν = 2., L1 = L1)
 sol0 = 1.1cos.(X) .* exp.(-0X.^2/(2*5^2))
 
 prob = BifurcationProblem(R_SH, sol0, parSH, (@optic _.λ); J = Jac_sp,
-    record_from_solution = (x, p) -> (n2 = norm(x), nw = normweighted(x), s = sum(x), s2 = x[end ÷ 2], s4 = x[end ÷ 4], s5 = x[end ÷ 5]),
+    record_from_solution = (x, p; k...) -> (n2 = norm(x), nw = normweighted(x), s = sum(x), s2 = x[end ÷ 2], s4 = x[end ÷ 4], s5 = x[end ÷ 5]),
     plot_solution = (x, p;kwargs...)->(plot!(X, x; ylabel="solution", label="", kwargs...))
     )
 ####################################################################################################
diff --git a/examples/TMModel.jl b/examples/TMModel.jl
index c57c6930..75e3ab60 100644
--- a/examples/TMModel.jl
+++ b/examples/TMModel.jl
@@ -18,14 +18,14 @@ end
 
 par_tm = (α = 1.5, τ = 0.013, J = 3.07, E0 = -2.0, τD = 0.200, U0 = 0.3, τF = 1.5, τS = 0.007)
 z0 = [0.238616, 0.982747, 0.367876 ]
-prob = BifurcationProblem(TMvf!, z0, par_tm, (@lens _.E0); record_from_solution = (x, p) -> (E = x[1], x = x[2], u = x[3]),)
+prob = BifurcationProblem(TMvf!, z0, par_tm, (@optic _.E0); record_from_solution = (x, p; k...) -> (E = x[1], x = x[2], u = x[3]),)
 
-opts_br = ContinuationPar(p_min = -10.0, p_max = -0., dsmax = 0.1, n_inversion = 8, nev = 3)
+opts_br = ContinuationPar(p_min = -10.0, p_max = 1., dsmax = 0.1, n_inversion = 8, nev = 3)
 br = continuation(prob, PALC(tangent = Bordered()), opts_br; plot = false, normC = norminf, bothside = true)
 
-BK.plot(br, plotfold=false)
+plot(br, plotfold=false)
 ####################################################################################################
-br_fold = BK.continuation(br, 1, (@optic _.α), ContinuationPar(br.contparams, p_min = 0.2, p_max = 5.), bothside = true)
+br_fold = BK.continuation(br, 2, (@optic _.α), ContinuationPar(br.contparams, p_min = 0.2, p_max = 5.), bothside = true)
 plot(br_fold)
 ####################################################################################################
 # continuation parameters
@@ -40,7 +40,7 @@ function plotSolution(x, p; k...)
     plot!(br; subplot = 1, putspecialptlegend = false)
 end
 
-args_po = (	record_from_solution = (x, p) -> begin
+args_po = (	record_from_solution = (x, p; k...) -> begin
         xtt = BK.get_periodic_orbit(p.prob, x, p.p)
         return (max = maximum(xtt[1,:]),
                 min = minimum(xtt[1,:]),
@@ -87,7 +87,7 @@ prob_ode = ODEProblem(TMvf!, copy(z0), (0., 1000.), par_tm; abstol = 1e-11, relt
 opts_po_cont = ContinuationPar(opts_br, ds= -0.0001, dsmin = 1e-4, max_steps = 120, newton_options = NewtonPar(tol = 1e-6, max_iterations = 7), tol_stability = 1e-8, detect_bifurcation = 2, plot_every_step = 10, save_sol_every_step=1)
 
 br_posh = @time continuation(
-    br, 4,
+    br, 5,
     # arguments for continuation
     opts_po_cont,
     # this is where we tell that we want Standard Shooting
@@ -105,7 +105,7 @@ plot(br_posh, br, markersize=3)
 opts_po_cont = ContinuationPar(opts_br, dsmax = 0.02, ds= 0.0001, max_steps = 50, newton_options = NewtonPar(tol = 1e-9, max_iterations=15), tol_stability = 1e-6, detect_bifurcation = 2, plot_every_step = 5)
 
 br_popsh = @time continuation(
-    br, 4,
+    br, 5,
     # arguments for continuation
     opts_po_cont,
     # this is where we tell that we want Poincaré Shooting
diff --git a/examples/brusselator.jl b/examples/brusselator.jl
index e63c5f11..d3f43ad9 100644
--- a/examples/brusselator.jl
+++ b/examples/brusselator.jl
@@ -96,7 +96,7 @@ prob = BifurcationProblem(Fbru!, sol0, par_bru, (@optic _.l);
         J = Jbru_sp, 
         # plot_solution = (x, p; kwargs...) -> plotsol(x; label="", kwargs... ), # for Plots.jl
         # plot_solution = (ax, x, p) -> plotsol(ax, x), # For Makie.jl
-        record_from_solution = (x, p) -> x[div(n,2)])
+        record_from_solution = (x, p; k...) -> x[div(n,2)])
 ####################################################################################################
 eigls = EigArpack(1.1, :LM)
 opts_br_eq = ContinuationPar(dsmin = 0.03, dsmax = 0.05, ds = 0.03, p_max = 1.9, detect_bifurcation = 3, nev = 21, plot_every_step = 50, newton_options = NewtonPar(eigsolver = eigls, tol = 1e-9), max_steps = 1060, n_inversion = 6, tol_bisection_eigenvalue = 1e-20, max_bisection_steps = 30)
diff --git a/examples/brusselatorShooting.jl b/examples/brusselatorShooting.jl
index 3ed0e1b3..6d6af2b2 100644
--- a/examples/brusselatorShooting.jl
+++ b/examples/brusselatorShooting.jl
@@ -90,7 +90,7 @@ sol0 = vcat(par_bru.α * ones(n), par_bru.β/par_bru.α * ones(n))
 probBif = BifurcationProblem(Fbru!, sol0, par_bru, (@optic _.l);
         J = Jbru_sp,
         plot_solution = (x, p; ax1 = 0, kwargs...) -> (plotsol(x; label="", kwargs... )),
-        record_from_solution = (x, p) -> x[div(n,2)]
+        record_from_solution = (x, p; k...) -> x[div(n,2)]
         )
 ####################################################################################################
 eigls = EigArpack(1.1, :LM)
@@ -177,7 +177,7 @@ br_po = @time continuation(deepcopy(probSh), outpo.u, PALC(),
         # (Base.display(contResult.eig[end].eigenvals) ;true),
     callback_newton = BK.cbMaxNorm(1.),
     plot_solution = (x, p; kwargs...) -> BK.plot_periodic_shooting!(x[1:end-1], length(1:dM:M); kwargs...),
-    record_from_solution = (u, p) -> u[end], normC = norminf)
+    record_from_solution = (u, p; k...) -> u[end], normC = norminf)
 
 ####################################################################################################
 # automatic branch switching with Shooting
diff --git a/examples/cGL2d.jl b/examples/cGL2d.jl
index 1da39e93..34e21c55 100644
--- a/examples/cGL2d.jl
+++ b/examples/cGL2d.jl
@@ -172,7 +172,7 @@ poTrap = PeriodicOrbitTrapProblem(re_make(prob, params = (@set par_cgl.r = r_hop
 
 ls0 = GMRESIterativeSolvers(N = 2n, reltol = 1e-9)#, Pl = lu(I + par_cgl.Δ))
 poTrapMF = setproperties(poTrap; linsolver = ls0)
-@set! poTrapMF.prob_vf.VF.J = (x, p) ->  (dx -> dFcgl(x, p, dx))
+@reset poTrapMF.prob_vf.VF.J = (x, p) ->  (dx -> dFcgl(x, p, dx))
 
 poTrap(orbitguess_f, @set par_cgl.r = r_hopf - 0.1) |> plot
 poTrapMF(orbitguess_f, @set par_cgl.r = r_hopf - 0.1) |> plot
@@ -205,7 +205,7 @@ ls = GMRESIterativeSolvers(verbose = false, reltol = 1e-3, N = size(Jpo,1), rest
 ls(Jpo, rand(ls.N))
 
 opt_po = @set opt_newton.verbose = true
-@set! opt_po.linsolver = ls
+@reset opt_po.linsolver = ls
 
 outpo_f = @time newton(poTrapMF, orbitguess_f, opt_po; normN = norminf)
 BK.converged(outpo_f) && printstyled(color=:red, "--> T = ", outpo_f.u[end], ", amplitude = ", BK.amplitude(outpo_f.u, Nx*Ny, M; ratio = 2),"\n")
@@ -219,7 +219,7 @@ br_po = @time continuation(poTrapMF, outpo_f.u, PALC(), opts_po_cont;
         verbosity = 3,
         plot = true,
         # plot_solution = (x, p;kwargs...) -> BK.plot_periodic_potrap(x, M, Nx, Ny; ratio = 2, kwargs...),
-        record_from_solution = (u, p) -> BK.getamplitude(poTrapMF, u, par_cgl; ratio = 2),
+        record_from_solution = (u, p; k...) -> BK.getamplitude(poTrapMF, u, par_cgl; ratio = 2),
         normC = norminf)
 
 branches = Any[br_pok2]
@@ -232,13 +232,15 @@ br_po = continuation(
     br, 1,
     # arguments for continuation
     opts_po_cont, poTrapMF;
-    ampfactor = 3.,
-    verbosity = 3, plot = true,
+    # ampfactor = 3.,
+    verbosity = 3, 
+    plot = true,
     # callback_newton = (x, f, J, res, iteration, itl, options; kwargs...) -> (println("--> amplitude = ", BK.amplitude(x, n, M; ratio = 2));true),
     finalise_solution = (z, tau, step, contResult; k...) ->
     (BK.haseigenvalues(contResult) && Base.display(contResult.eig[end].eigenvals) ;true),
     plot_solution = (x, p; kwargs...) -> BK.plot_periodic_potrap(x, M, Nx, Ny; ratio = 2, kwargs...),
-    record_from_solution = (u, p) -> BK.amplitude(u, Nx*Ny, M; ratio = 2), normC = norminf)
+    record_from_solution = (u, p; k...) -> BK.amplitude(u, Nx*Ny, M; ratio = 2), 
+    normC = norminf)
 ####################################################################################################
 # Experimental, full Inplace
 @views function NL!(f, u, p, t = 0.)
@@ -302,7 +304,7 @@ out_ = similar(sol0f)
 @time Fcgl!(out_, sol0f, par_cgl)
 @time dFcgl!(out_, sol0f, par_cgl, sol0f)
 
-probInp = BifurcationProblem(Fcgl!, vec(sol0), (@set par_cgl.r = r_hopf - 0.01), (@lens _.r); J = dFcgl!, inplace = true)
+probInp = BifurcationProblem(Fcgl!, vec(sol0), (@set par_cgl.r = r_hopf - 0.01), (@optic _.r); J = dFcgl!, inplace = true)
 
 ls = GMRESIterativeSolvers(verbose = false, reltol = 1e-3, N = size(Jpo,1), restart = 40, maxiter = 50, Pl = Precilu, log=true)
 ls(Jpo, rand(ls.N))
diff --git a/examples/carrier.jl b/examples/carrier.jl
index 5eb56ca3..5303f9b5 100644
--- a/examples/carrier.jl
+++ b/examples/carrier.jl
@@ -36,7 +36,7 @@ X = LinRange(-1,1,N)
 dx = X[2] - X[1]
 par_car = (ϵ = 0.7, X = X, dx = dx)
 sol = -(1 .- par_car.X.^2)
-recordFromSolution(x, p) = (x[2]-x[1]) * sum(x->x^2, x)
+recordFromSolution(x, p; k...) = (x[2]-x[1]) * sum(x->x^2, x)
 
 prob = BK.BifurcationProblem(F_carr, zeros(N), par_car, (@optic _.ϵ);
     J = Jac_carr,
diff --git a/examples/chan-af.jl b/examples/chan-af.jl
index 87cbc3d4..b1004718 100644
--- a/examples/chan-af.jl
+++ b/examples/chan-af.jl
@@ -59,7 +59,9 @@ end
 sol0 = Fun( x -> x * (1-x), Interval(0.0, 1.0))
 const Δ = Derivative(sol0.space, 2);
 par_af = (α = 3., β = 0.01)
-prob = BifurcationProblem(F_chan, sol0, par_af, (@optic _.α); J = Jac_chan, plot_solution = (x, p; kwargs...) -> plot!(x; label = "l = $(length(x))", kwargs...))
+prob = BifurcationProblem(F_chan, sol0, par_af, (@optic _.α); 
+            J = Jac_chan, 
+            plot_solution = (x, p; kwargs...) -> plot!(x; label = "l = $(length(x))", kwargs...))
 
 optnew = NewtonPar(tol = 1e-12, verbose = true)
 sol = @time BK.newton(prob, optnew, normN = norminf)
diff --git a/examples/codim2PO-sh-mf.jl b/examples/codim2PO-sh-mf.jl
index 852fd797..a03a84c1 100644
--- a/examples/codim2PO-sh-mf.jl
+++ b/examples/codim2PO-sh-mf.jl
@@ -21,10 +21,10 @@ par_pop = ComponentArray( K = 1., r = 2π, a = 4π, b0 = 0.25, e = 1., d = 2π,
 
 z0 = [0.1,0.1,1,0]
 
-prob = BifurcationProblem(Pop!, z0, par_pop, (@optic _.b0); record_from_solution = (x, p) -> (x = x[1], y = x[2], u = x[3]))
+prob = BifurcationProblem(Pop!, z0, par_pop, (@optic _.b0); record_from_solution = (x, p; k...) -> (x = x[1], y = x[2], u = x[3]))
 
 opts_br = ContinuationPar(p_min = 0., p_max = 20.0, ds = 0.002, dsmax = 0.01, n_inversion = 6, detect_bifurcation = 3, max_bisection_steps = 25, nev = 4, max_steps = 20000)
-@set! opts_br.newton_options.verbose = true
+@reset opts_br.newton_options.verbose = true
 
 ################################################################################
 using DifferentialEquations
@@ -37,7 +37,7 @@ sol = solve(prob_de, Rodas5())
 
 plot(sol)
 ################################################################################
-argspo = (record_from_solution = (x, p) -> begin
+argspo = (record_from_solution = (x, p; k...) -> begin
         xtt = BK.get_periodic_orbit(p.prob, x, set(getparams(p.prob), BK.getlens(p.prob), p.p))
         return (max = maximum(xtt[1,:]),
                 min = minimum(xtt[1,:]),
@@ -92,7 +92,7 @@ res1 = AD.pullback_function(AD.ZygoteBackend(), x->flow(x, prob_de,1), sol0_f)(s
 
 AD.pullback_function(AD.FiniteDifferencesBackend(), z -> probsh(z, getparams(probsh)), cish)(cish)[1]
 
-@set! probsh.flow.vjp = (x,p,dx,tm) -> AD.pullback_function(AD.ZygoteBackend(), z->flow(z, prob_de,tm,p), x)(dx)[1]
+@reset probsh.flow.vjp = (x,p,dx,tm) -> AD.pullback_function(AD.ZygoteBackend(), z->flow(z, prob_de,tm,p), x)(dx)[1]
 
 lspo = GMRESIterativeSolvers(verbose = false, N = length(cish), abstol = 1e-12, reltol = 1e-10)
 eigpo = EigKrylovKit(x₀ = rand(4))
@@ -103,7 +103,7 @@ _sol = BK.get_periodic_orbit(probsh, solpo.u, sol.prob.p)
 plot(_sol.t, _sol[1:2,:]')
 
 opts_po_cont = setproperties(opts_br, max_steps = 50, save_eigenvectors = true, detect_loop = true, tol_stability = 1e-3, newton_options = optnpo)
-@set! opts_po_cont.newton_options.verbose = true
+@reset opts_po_cont.newton_options.verbose = true
 br_fold_sh = continuation(probsh, cish, PALC(tangent = Bordered()), opts_po_cont;
     verbosity = 3, plot = true,
     linear_algo = MatrixFreeBLS(@set lspo.N = lspo.N+1),
@@ -121,10 +121,10 @@ brpo_pd_sh = continuation(probsh2, cish, PALC(), opts_po_cont;
 
 # codim 2 Fold
 opts_posh_fold = ContinuationPar(br_fold_sh.contparams, detect_bifurcation = 0, max_steps = 20, p_min = 0.01, p_max = 1.2)
-@set! opts_posh_fold.newton_options.tol = 1e-9
+@reset opts_posh_fold.newton_options.tol = 1e-9
 
 # use this option for jacobian_ma = :finiteDifferencesMF, otherwise do not
-@set! opts_posh_fold.newton_options.linsolver.solver.N = opts_posh_fold.newton_options.linsolver.solver.N+1
+@reset opts_posh_fold.newton_options.linsolver.solver.N = opts_posh_fold.newton_options.linsolver.solver.N+1
 fold_po_sh1 = continuation(br_fold_sh, 2, (@optic _.ϵ), opts_posh_fold;
     verbosity = 2, plot = true,
     detect_codim2_bifurcation = 0,
@@ -157,9 +157,9 @@ plot(fold_po_sh1, fold_po_sh2, branchlabel = ["FOLD", "FOLD"])
 
 # codim 2 PD
 opts_posh_pd = ContinuationPar(brpo_pd_sh.contparams, detect_bifurcation = 3, max_steps = 40, p_min = -1.)
-@set! opts_posh_pd.newton_options.tol = 1e-8
+@reset opts_posh_pd.newton_options.tol = 1e-8
 # use this option for jacobian_ma = :finiteDifferencesMF, otherwise do not
-# @set! opts_posh_pd.newton_options.linsolver.solver.N = opts_posh_pd.newton_options.linsolver.solver.N+1
+# @reset opts_posh_pd.newton_options.linsolver.solver.N = opts_posh_pd.newton_options.linsolver.solver.N+1
 pd_po_sh = continuation(brpo_pd_sh, 1, (@optic _.b0), opts_posh_pd;
     verbosity = 2, plot = true,
     detect_codim2_bifurcation = 0,
@@ -190,7 +190,7 @@ probshns, ci = generate_ci_problem( ShootingProblem(M=3), re_make(prob, params =
             # jacobian = BK.FiniteDifferencesMF()
             )
 
-@set! probshns.flow.vjp = (x,p,dx,tm) -> AD.pullback_function(AD.ZygoteBackend(), z->flow(z, prob_de,tm,p), x)(dx)[1]
+@reset probshns.flow.vjp = (x,p,dx,tm) -> AD.pullback_function(AD.ZygoteBackend(), z->flow(z, prob_de,tm,p), x)(dx)[1]
 
 brpo_ns = continuation(probshns, ci, PALC(), ContinuationPar(opts_po_cont; max_steps = 50, ds = -0.001);
     verbosity = 3, plot = true,
@@ -205,7 +205,7 @@ ns = get_normal_form(brpo_ns, 1)
 # codim 2 NS
 using AbbreviatedStackTraces
 opts_posh_ns = ContinuationPar(brpo_ns.contparams, detect_bifurcation = 0, max_steps = 100, p_min = -0., p_max = 1.2)
-@set! opts_posh_ns.newton_options.tol = 1e-9
+@reset opts_posh_ns.newton_options.tol = 1e-9
 ns_po_sh = continuation(brpo_ns, 1, (@optic _.ϵ), opts_posh_ns;
         verbosity = 2, plot = true,
         detect_codim2_bifurcation = 0,
diff --git a/examples/codim2PO.jl b/examples/codim2PO.jl
index 9e1fa3a4..3b47ad35 100644
--- a/examples/codim2PO.jl
+++ b/examples/codim2PO.jl
@@ -36,7 +36,7 @@ sol = solve(prob_de, Rodas5())
 
 plot(sol)
 ################################################################################
-function recordFromSolution(x, p)
+function recordFromSolution(x, p; k...)
     xtt = BK.get_periodic_orbit(p.prob, x, p.p)
     _max = maximum(xtt[1,:])
     _min = minimum(xtt[1,:])
@@ -148,19 +148,19 @@ _sol = BK.get_periodic_orbit(probcoll, solpo.u,1)
 plot(_sol.t, _sol[1:2,:]')
 
 opts_po_cont = ContinuationPar(opts_br, max_steps = 50, tol_stability = 1e-8, n_inversion = 6)
-@set! opts_po_cont.newton_options.verbose = true
+@reset opts_po_cont.newton_options.verbose = true
 brpo_fold = continuation(probcoll, ci, PALC(), opts_po_cont;
     verbosity = 3, plot = true,
     argspo...
     )
-pd = get_normal_form(brpo_fold, 1; prm = true)
+pd = get_normal_form(brpo_fold, 1; prm = false)
 
 prob2 = @set probcoll.prob_vf.lens = @optic _.ϵ
 brpo_pd = continuation(prob2, ci, PALC(), ContinuationPar(opts_po_cont, dsmax = 5e-3, max_steps=250);
     verbosity = 3, plot = true,
     argspo...
     )
-pt = get_normal_form(brpo_pd, 1, prm = true)
+pt = get_normal_form(brpo_pd, 1, prm = false)
 #########
 # pd branch switching
 brpo_pd_sw = continuation(deepcopy(brpo_pd), 1,
@@ -180,8 +180,8 @@ opts_pocoll_fold = ContinuationPar(brpo_fold.contparams,
                                   p_max = 1.2,
                                   n_inversion = 4,
                                   plot_every_step = 10)
-@set! opts_pocoll_fold.newton_options.tol = 1e-12
-fold_po_coll1 = continuation(brpo_fold, 1, (@lens _.ϵ), opts_pocoll_fold;
+@reset opts_pocoll_fold.newton_options.tol = 1e-12
+fold_po_coll1 = continuation(brpo_fold, 1, (@optic _.ϵ), opts_pocoll_fold;
         verbosity = 3, plot = true,
         detect_codim2_bifurcation = 0,
         start_with_eigen = false,
diff --git a/examples/mittleman.jl b/examples/mittleman.jl
index 60560e92..ab551140 100644
--- a/examples/mittleman.jl
+++ b/examples/mittleman.jl
@@ -75,7 +75,7 @@ w .-= minimum(w)
 
 prob = BK.BifurcationProblem(Fmit!, sol0, par_mit, (@optic _.λ);
         J = JFmit,
-        record_from_solution = (x, p) -> (nw = normbratu(x), n2 = norm(x), n∞ = norminf(x)),
+        record_from_solution = (x, p; k...) -> (nw = normbratu(x), n2 = norm(x), n∞ = norminf(x)),
         plot_solution = (x, p; k...) -> plotsol!(x ; k...),
         # plot_solution = (ax, x, p; ax1 = nothing, k...) -> plotsol!(ax, x ; k...),
         issymmetric = true)
@@ -118,16 +118,16 @@ opts_br = ContinuationPar(dsmin = 0.0001, dsmax = 0.04, ds = 0.005, p_max = 3.5,
 
 br = @time continuation(prob, PALC(), opts_br; kwargsC...)
 
-BK.plot(br)
+plot(br)
 ####################################################################################################
 # automatic branch switching
 br1 = continuation(br, 3; kwargsC...)
 
-BK.plot(br, br1, plotfold=false)
+plot(br, br1, plotfold=false)
 
 br2 = continuation(br1, 1; kwargsC...)
 
-BK.plot(br, br1, br2, plotfold=false)
+plot(br, br1, br2, plotfold=false)
 ####################################################################################################
 # bifurcation diagram
 function optionsCont(x,p,l; opt = opts_br)
diff --git a/examples/pd-1d.jl b/examples/pd-1d.jl
index edece36e..489e584a 100644
--- a/examples/pd-1d.jl
+++ b/examples/pd-1d.jl
@@ -61,7 +61,7 @@ u0 = cos.(2X)
 solc0 = vcat(u0, u0)
 
 probBif = BifurcationProblem(Fbr, solc0, par_br, (@optic _.C) ;J = Jbr,
-    record_from_solution = (x, p) -> norm(x, Inf),
+    record_from_solution = (x, p; k...) -> norm(x, Inf),
     plot_solution = (x, p; kwargs...) -> plot!(x[1:end÷2];label="",ylabel ="u", kwargs...))
 ####################################################################################################
 # eigls = DefaultEig()
@@ -101,7 +101,7 @@ br_po = @time continuation(
     plot = true,
     callback_newton = BK.cbMaxNorm(1e2),
     plot_solution = (x, p;kwargs...) ->  (heatmap!(reshape(x[1:end-1], 2*N, M)'; ylabel="time", color=:viridis, kwargs...);plot!(br, subplot=1)),
-    record_from_solution = (u, p) -> (max = maximum(u[1:end-1]), period = u[end]),#BK.maximumPOTrap(u, N, M; ratio = 2),
+    record_from_solution = (u, p; k...) -> (max = maximum(u[1:end-1]), period = u[end]),#BK.maximumPOTrap(u, N, M; ratio = 2),
     normC = norminf)
 
 plot(br, br_po, label = "")
@@ -121,7 +121,7 @@ br_po_pd = @time continuation(
     # jacobianPO = :FullSparseInplace,
     # jacobianPO = :BorderedSparseInplace,
     plot_solution = (x, p;kwargs...) ->  (heatmap!(reshape(x[1:end-1], 2*N, M)'; ylabel="time", color=:viridis, kwargs...);plot!(br_po, subplot=1)),
-    record_from_solution = (u, p) -> (max = maximum(u[1:end-1]), period = u[end]),#BK.maximumPOTrap(u, N, M; ratio = 2),
+    record_from_solution = (u, p; k...) -> (max = maximum(u[1:end-1]), period = u[end]),#BK.maximumPOTrap(u, N, M; ratio = 2),
     normC = norminf)
 
 plot(br, br_po, br_po_pd, label = "")
@@ -167,7 +167,7 @@ br_po_sh = @time continuation(probSh, outposh.u, PALC(), optcontpo;
     finalise_solution = (z, tau, step, contResult; kw...) ->
         (BK.haseigenvalues(contResult) && Base.display(contResult.eig[end].eigenvals) ;true),
     plot_solution = (x, p; kwargs...) -> BK.plot_periodic_shooting!(x[1:end-1], 1; kwargs...),
-    record_from_solution = (u, p) -> BK.getmaximum(probSh, u, (@set par_br_hopf.C = p.p); ratio = 2), normC = norminf)
+    record_from_solution = (u, p; k...) -> BK.getmaximum(probSh, u, (@set par_br_hopf.C = p.p); ratio = 2), normC = norminf)
 
 # branches = [br_po_sh]
 # push!(branches, br_po_sh)
@@ -251,7 +251,7 @@ br_po = @time continuation(
     finalise_solution = (z, tau, step, contResult; kw...) ->
         (BK.haseigenvalues(contResult) && Base.display(contResult.eig[end].eigenvals) ;true),
     plot_solution = (x, p; kwargs...) -> (BK.plot_periodic_shooting!(x[1:end-1], 1; kwargs...);plot!(br, subplot=1)),
-    record_from_solution = (u, p) -> BK.getmaximum(probPO, u, (@set par_br.C = p.p); ratio = 2),
+    record_from_solution = (u, p; k...) -> BK.getmaximum(probPO, u, (@set par_br.C = p.p); ratio = 2),
     normC = norminf)
 
 plot(br, br_po, label = "")
@@ -273,7 +273,7 @@ br_po_pd = BK.continuation(br_po, 1, setproperties(br_po.contparams, detect_bifu
         heatmap!(outt[:,:]'; color = :viridis, subplot = 3)
         plot!(br_po; legend=false, subplot=1)
     end,
-    record_from_solution = (u, p) -> (BK.getmaximum(p.prob, u, (@set par_br_hopf.C = p.p); ratio = 2)), normC = norminf
+    record_from_solution = (u, p; k...) -> (BK.getmaximum(p.prob, u, (@set par_br_hopf.C = p.p); ratio = 2)), normC = norminf
     )
 
 plot(br_po, br_po_pd, legend=false)
@@ -296,7 +296,7 @@ br_po_pdcodim2 = @time continuation(
     finalise_solution = (z, tau, step, contResult; kw...) ->
         (BK.haseigenvalues(contResult) && Base.display(contResult.eig[end].eigenvals) ;true),
     plot_solution = (x, p; kwargs...) -> (BK.plot_periodic_shooting!(x[1:end-1], 1; kwargs...);plot!(br, subplot=1)),
-    record_from_solution = (u, p) -> BK.getmaximum(probPO, u, (@set par_br.C = p.p); ratio = 2),
+    record_from_solution = (u, p; k...) -> BK.getmaximum(probPO, u, (@set par_br.C = p.p); ratio = 2),
     normC = norminf)
 ####################################################################################################
 # aBS Poincare Shooting
diff --git a/src/Continuation.jl b/src/Continuation.jl
index b7a22740..c0eba489 100644
--- a/src/Continuation.jl
+++ b/src/Continuation.jl
@@ -245,7 +245,7 @@ end
 function get_state_summary(it, state::ContState{Tv, T, Teigvals}) where {Tv, T, Teigvals}
     x = getx(state)
     p = getp(state)
-    pt = record_from_solution(it)(x, p)
+    pt = record_from_solution(it)(x, p; iter = it, state = state)
     stable = Teigvals!=Nothing ? is_stable(state) : nothing
     return mergefromuser(pt, 
                         (param = p,
@@ -304,7 +304,7 @@ function ContResult(it::AbstractContinuationIterable,
                     state::AbstractContinuationState)
     x0 = _copy(getx(state))
     p0 = getp(state)
-    pt = record_from_solution(it)(x0, p0)
+    pt = record_from_solution(it)(x0, p0; iter = it, state = state)
     return _contresult(it, state,
                         pt,
                         get_state_summary(it, state), 
diff --git a/src/Problems.jl b/src/Problems.jl
index eef1851f..68f866de 100644
--- a/src/Problems.jl
+++ b/src/Problems.jl
@@ -16,7 +16,7 @@ Determine if the vector field is of the form `f!(out,z,p)`.
 """
 function _isinplace(f)
     m = minimum(numargs(f))
-    @assert 1 < m < 4 "You have too many/few arguments in your vector field F. It should be of the form `F(x,p)` or `F!(x,p)`."
+    @assert 1 < m < 4 "You have too many/few arguments in your vector field F. It should be of the form `F(x, p)` or `F!(x, p)`."
     return m == 3
 end
 
@@ -140,7 +140,7 @@ for (op, at) in (
                 lens::Tl
                 "user function to plot solutions during continuation. Signature: `plot_solution(x, p; kwargs...)` for Plot.jl and `plot_solution(ax, x, p; kwargs...)` for the Makie package(s)."
                 plotSolution::Tplot
-                "`record_from_solution = (x, p) -> norm(x)` function used record a few indicators about the solution. It could be `norm` or `(x, p) -> x[1]`. This is also useful when saving several huge vectors is not possible for memory reasons (for example on GPU). This function can return pretty much everything but you should keep it small. For example, you can do `(x, p) -> (x1 = x[1], x2 = x[2], nrm = norm(x))` or simply `(x, p) -> (sum(x), 1)`. This will be stored in `contres.branch` where `contres::ContResult` is the continuation curve of the bifurcation problem. Finally, the first component is used for plotting in the continuation curve."
+                "`record_from_solution = (x, p; k...) -> norm(x)` function used record a few indicators about the solution. It could be `norm` or `(x, p) -> x[1]`. This is also useful when saving several huge vectors is not possible for memory reasons (for example on GPU). This function can return pretty much everything but you should keep it small. For example, you can do `(x, p; k...) -> (x1 = x[1], x2 = x[2], nrm = norm(x))` or simply `(x, p; k...) -> (sum(x), 1)`. This will be stored in `contres.branch` where `contres::ContResult` is the continuation curve of the bifurcation problem. Finally, the first component is used for plotting in the continuation curve."
                 recordFromSolution::Trec
                 "function to save the full solution on the branch. Some problem are mutable (like periodic orbit functional with adaptive mesh) and this function allows to save the state of the problem along with the solution itself. Signature `save_solution(x, p)`"
                 save_solution::Tgets
diff --git a/src/periodicorbit/BifurcationPoints.jl b/src/periodicorbit/BifurcationPoints.jl
index 3eb92df1..41b7d4ae 100644
--- a/src/periodicorbit/BifurcationPoints.jl
+++ b/src/periodicorbit/BifurcationPoints.jl
@@ -131,7 +131,7 @@ function Base.show(io::IO, ns::NeimarkSackerPO)
     println(io, "├─ Period at the periodic orbit T ≈ ", abs(ns.T))
     println(io, "├─ Second period of the bifurcated torus ≈ ", abs(2pi*ns.ω*ns.T))
     if ns.prm
-        println(io, "├─ Normal form z -> z⋅eⁱᶿ(1 + a⋅δp + b⋅|z|²)")
+        println(io, "├─ Normal form z ─▶ z⋅eⁱᶿ(1 + a⋅δp + b⋅|z|²)")
     else
         println(io, "├─ Normal form:\n├\t∂τ = 1 + a⋅|ξ|²\n├\t∂ξ = iθ/T⋅ξ + d⋅ξ⋅|ξ|²")
     end
diff --git a/src/periodicorbit/PeriodicOrbitTrapeze.jl b/src/periodicorbit/PeriodicOrbitTrapeze.jl
index 09b0ff90..b3161943 100644
--- a/src/periodicorbit/PeriodicOrbitTrapeze.jl
+++ b/src/periodicorbit/PeriodicOrbitTrapeze.jl
@@ -1069,7 +1069,7 @@ function continuation(prob::PeriodicOrbitTrapProblem,
                     orbitguess,
                     alg::AbstractContinuationAlgorithm,
                     _contParams::ContinuationPar;
-                    record_from_solution = (u, p) -> (period = u[end],),
+                    record_from_solution = (u, p; k...) -> (period = u[end],),
                     linear_algo = nothing,
                     kwargs...)
     _linear_algo = isnothing(linear_algo) ?  BorderingBLS(solver = _contParams.newton_options.linsolver, check_precision = false) : linear_algo
diff --git a/test/COModel.jl b/test/COModel.jl
index 0ec594d7..436d511b 100644
--- a/test/COModel.jl
+++ b/test/COModel.jl
@@ -20,7 +20,7 @@ par_com = (q1 = 2.5, q2 = 1., q3 = 10., q4 = 0.0675, q5 = 1., q6 = 0.1, k = 0.4)
 z0 = [0.001137, 0.891483, 0.062345]
 
 prob = BifurcationProblem(COm, z0, par_com, (@optic _.q2);
-        record_from_solution = (x, p) -> (x = x[1], y = x[2], s = x[3]))
+        record_from_solution = (x, p; k...) -> (x = x[1], y = x[2], s = x[3]))
 
 opts_br = ContinuationPar(p_min = 0.5, p_max = 2.3, ds = 0.002, dsmax = 0.01, n_inversion = 6, detect_bifurcation = 3, max_bisection_steps = 25, nev = 3, max_steps = 100)
 
diff --git a/test/codim2.jl b/test/codim2.jl
index 08f3803e..8c204bd5 100644
--- a/test/codim2.jl
+++ b/test/codim2.jl
@@ -18,7 +18,7 @@ end
 par_com = (q1 = 2.5, q2 = 2.0, q3 = 10., q4 = 0.0675, q5 = 1., q6 = 0.1, k = 0.4)
 z0 = [0.07,0.2,05]
 
-prob = BifurcationProblem(COm, z0, par_com, (@optic _.q2); record_from_solution = (x, p) -> (x = x[1], y = x[2], s = x[3]))
+prob = BifurcationProblem(COm, z0, par_com, (@optic _.q2); record_from_solution = (x, p; k...) -> (x = x[1], y = x[2], s = x[3]))
 
 opts_br = ContinuationPar(p_min = 0.6, p_max = 2.5, ds = 0.002, dsmax = 0.01, n_inversion = 4, detect_bifurcation = 3, max_bisection_steps = 25, nev = 2, max_steps = 20000)
 
diff --git a/test/codim2PO-OColl.jl b/test/codim2PO-OColl.jl
index 44d3c2d9..9d1aab30 100644
--- a/test/codim2PO-OColl.jl
+++ b/test/codim2PO-OColl.jl
@@ -19,7 +19,7 @@ par_pop = ( K = 1., r = 2π, a = 4π, b0 = 0.25, e = 1., d = 2π, ϵ = 0.2, )
 
 z0 = [0.1,0.1,1,0]
 
-prob = BifurcationProblem(Pop!, z0, par_pop, (@optic _.b0); record_from_solution = (x, p) -> (x = x[1], y = x[2], u = x[3]))
+prob = BifurcationProblem(Pop!, z0, par_pop, (@optic _.b0); record_from_solution = (x, p; k...) -> (x = x[1], y = x[2], u = x[3]))
 
 opts_br = ContinuationPar(p_min = 0., p_max = 20.0, ds = 0.002, dsmax = 0.01, n_inversion = 6, detect_bifurcation = 3, max_bisection_steps = 25, nev = 4, max_steps = 20000)
 @reset opts_br.newton_options.verbose = false
@@ -32,7 +32,7 @@ sol = solve(prob_de, alg)
 prob_de = ODEProblem(Pop!, sol.u[end], (0,5.), par_pop, reltol = 1e-8, abstol = 1e-10)
 sol = solve(prob_de, Rodas5())
 ################################################################################
-argspo = (record_from_solution = (x, p) -> begin
+argspo = (record_from_solution = (x, p; k...) -> begin
         xtt = BK.get_periodic_orbit(p.prob, x, p.p)
         return (max = maximum(xtt[1,:]),
                 min = minimum(xtt[1,:]),
diff --git a/test/codim2PO-shooting-mf.jl b/test/codim2PO-shooting-mf.jl
index 8c1380a8..622a5ba8 100644
--- a/test/codim2PO-shooting-mf.jl
+++ b/test/codim2PO-shooting-mf.jl
@@ -36,7 +36,7 @@ prob_de = ODEProblem(Pop!, sol.u[end], (0,5.), par_pop, reltol = 1e-10, abstol =
 sol = solve(prob_de, Rodas5())
 ################################################################################
 @info "plotting function"
-argspo = (record_from_solution = (x, p) -> begin
+argspo = (record_from_solution = (x, p; k...) -> begin
         xtt = BK.get_periodic_orbit(p.prob, x, p.p)
         return (max = maximum(xtt[1,:]),
                 min = minimum(xtt[1,:]),
diff --git a/test/codim2PO-shooting.jl b/test/codim2PO-shooting.jl
index ca0eee5d..22b86d64 100644
--- a/test/codim2PO-shooting.jl
+++ b/test/codim2PO-shooting.jl
@@ -19,7 +19,7 @@ par_pop = ( K = 1., r = 2π, a = 4π, b0 = 0.25, e = 1., d = 2π, ϵ = 0.2, )
 
 z0 = [0.1,0.1,1,0]
 
-prob = BifurcationProblem(Pop!, z0, par_pop, (@optic _.b0); record_from_solution = (x, p) -> (x = x[1], y = x[2], u = x[3]))
+prob = BifurcationProblem(Pop!, z0, par_pop, (@optic _.b0); record_from_solution = (x, p; k...) -> (x = x[1], y = x[2], u = x[3]))
 
 opts_br = ContinuationPar(p_min = 0., p_max = 20.0, ds = 0.002, dsmax = 0.01, n_inversion = 6, detect_bifurcation = 3, max_bisection_steps = 25, nev = 4, max_steps = 200)
 ################################################################################
@@ -30,7 +30,7 @@ sol = solve(prob_de, alg)
 prob_de = ODEProblem(Pop!, sol.u[end], (0,5.), par_pop, reltol = 1e-8, abstol = 1e-10)
 sol = solve(prob_de, Rodas5())
 ################################################################################
-argspo = (record_from_solution = (x, p) -> begin
+argspo = (record_from_solution = (x, p; k...) -> begin
         xtt = BK.get_periodic_orbit(p.prob, x, p.p)
         return (max = maximum(xtt[1,:]),
                 min = minimum(xtt[1,:]),
diff --git a/test/event.jl b/test/event.jl
index 765839c7..a95ca0fe 100644
--- a/test/event.jl
+++ b/test/event.jl
@@ -72,7 +72,7 @@ par = (p1 = -3., p2=-3., k=3)
 opts0 = ContinuationPar(dsmax = 0.1, ds = 0.001, max_steps = 1000, p_min = -3., p_max = 4.0, save_sol_every_step = 1, newton_options = NewtonPar(tol = 1e-10, verbose = false, max_iterations = 5), detect_bifurcation = 3, detect_event = 0, n_inversion = 8, dsmin_bisection = 1e-9, max_bisection_steps = 15, detect_fold=false, plot_every_step = 10)
 
 prob = BK.BifurcationProblem(Feve, -2ones(2), par, (@optic _.p1);
-        record_from_solution = (x, p) -> x[1])
+        record_from_solution = (x, p; k...) -> x[1])
 
 br0 = continuation(prob, PALC(), opts0;
     plot = false, verbosity = 0,
@@ -90,7 +90,7 @@ br = continuation(prob, PALC(), opts;
 # pretty_table(br.branch[1:40])
 
 # arguments for continuation
-args = (BK.re_make(prob; record_from_solution = (x,p) -> x[1]), PALC(), opts)
+args = (BK.re_make(prob; record_from_solution = (x,p; k...) -> x[1]), PALC(), opts)
 kwargs = (plot = false, verbosity = 0, linear_algo = MatrixBLS(),)
 
 br = continuation(args...; kwargs...,
diff --git a/test/lorenz84.jl b/test/lorenz84.jl
index 21072917..eb614873 100644
--- a/test/lorenz84.jl
+++ b/test/lorenz84.jl
@@ -29,8 +29,8 @@ opts_br = ContinuationPar(p_min = -1.5, p_max = 3.0, ds = 0.001, dsmax = 0.025,
 
 z0 =  [2.9787004394953343, -0.03868302503393752,  0.058232737694740085, -0.02105288273117459]
 
-recordFromSolutionLor(u::AbstractVector, p) = (X = u[1], Y = u[2], Z = u[3], U = u[4])
-recordFromSolutionLor(u::BorderedArray, p) = recordFromSolutionLor(u.u, p)
+recordFromSolutionLor(u::AbstractVector, p; k...) = (X = u[1], Y = u[2], Z = u[3], U = u[4])
+recordFromSolutionLor(u::BorderedArray, p; k...) = recordFromSolutionLor(u.u, p)
 
 prob = BK.BifurcationProblem(Lor, z0, parlor, (@optic _.F);
     record_from_solution = recordFromSolutionLor,)
diff --git a/test/simple_continuation.jl b/test/simple_continuation.jl
index 31c0322f..0146a379 100644
--- a/test/simple_continuation.jl
+++ b/test/simple_continuation.jl
@@ -168,7 +168,7 @@ BK.hassolution(Branch(br1, nothing)) # test the method
 br1.param
 BK.getparams(br1)
 
-@reset prob.recordFromSolution = (x,p) -> norm(x,2)
+@reset prob.recordFromSolution = (x,p; k...) -> norm(x,2)
 br2 = continuation(prob, PALC(), opts)
 
 # test for different norms
@@ -193,7 +193,7 @@ br5a = continuation(prob, PALC(), opts, finalise_solution = finalise_solution)
 br6 = continuation(prob, PALC(tangent = Secant()), opts)
 
 optsnat = setproperties(opts; ds = 0.001, dsmax = 0.1, dsmin = 0.0001)
-br7 = continuation((@set prob.recordFromSolution = (x,p)->x[1]), Natural(), optsnat)
+br7 = continuation((@set prob.recordFromSolution = (x,p;k...)->x[1]), Natural(), optsnat)
 
 # tangent prediction with Bordered predictor
 br8 = continuation(prob, PALC(tangent = Bordered()), opts)
diff --git a/test/stuartLandauTrap.jl b/test/stuartLandauTrap.jl
index 70929555..e1deab38 100644
--- a/test/stuartLandauTrap.jl
+++ b/test/stuartLandauTrap.jl
@@ -65,7 +65,7 @@ for (ind, jacobianPO) in enumerate((:Dense, :DenseAD, :FullLU, :BorderedLU, :Ful
         PALC(),    (@set opts_po_cont.newton_options.linsolver = _ls);
         verbosity = 0, plot = false,
         linear_algo = BorderingBLS(solver = _ls, check_precision = false),
-        record_from_solution = (u, p) -> BK.getamplitude(poTrap, u, par_hopf; ratio = 1), normC = norminf)
+        record_from_solution = (u, p; k...) -> BK.getamplitude(poTrap, u, par_hopf; ratio = 1), normC = norminf)
 
     BK.get_periodic_orbit(br_po, 1)
 end
diff --git a/test/test-cont-non-vector.jl b/test/test-cont-non-vector.jl
index de977eec..fc6b0dfd 100644
--- a/test/test-cont-non-vector.jl
+++ b/test/test-cont-non-vector.jl
@@ -22,7 +22,7 @@ end
 
 opt_newton0 = NewtonPar(tol = 1e-11, verbose = false)
 prob = BK.BifurcationProblem(F0, [0.8], 1., (@optic _);
-        record_from_solution = (x, p) -> x[1],
+        record_from_solution = (x, p;k...) -> x[1],
         J = (x, r) -> diagm(0 => 1 .- 3 .* x.^2),
         Jᵗ = (x, r) -> diagm(0 => 1 .- 3 .* x.^2),
         d2F = (x, r, v1, v2) -> -6 .* x .* v1 .* v2,)
@@ -72,7 +72,7 @@ end
 sol0 = BorderedArray([0.8], 0.0)
 
 opt_newton = NewtonPar(tol = 1e-11, verbose = false, linsolver = linsolveBd())
-prob = BK.BifurcationProblem(Fb, sol0, (1., 1.), (@optic _[1]); J = (x, r) -> Jacobian(x, r[1], r[2]), record_from_solution = (x,p) -> x.u[1])
+prob = BK.BifurcationProblem(Fb, sol0, (1., 1.), (@optic _[1]); J = (x, r) -> Jacobian(x, r[1], r[2]), record_from_solution = (x,p;k...) -> x.u[1])
 sol = newton(prob, opt_newton)
 @test BK.converged(sol)
 
diff --git a/test/testHopfMA.jl b/test/testHopfMA.jl
index 9f9724db..4eda897e 100644
--- a/test/testHopfMA.jl
+++ b/test/testHopfMA.jl
@@ -71,7 +71,7 @@ par_bru = (α = 2., β = 5.45, D1 = 0.008, D2 = 0.004, l = 0.3)
 sol0 = vcat(par_bru.α * ones(n), par_bru.β/par_bru.α * ones(n))
 prob = BifurcationKit.BifurcationProblem(Fbru, sol0, par_bru, (@optic _.l);
         J = Jbru_ana,
-        record_from_solution = (x, p) -> norminf(x))
+        record_from_solution = (x, p; k...) -> norminf(x))
 
 # test that the jacobian is well computed
 @test Jbru_sp(sol0, par_bru) - Jbru_ana(sol0, par_bru) |> sparse |> nnz == 0
diff --git a/test/testJacobianFoldDeflation.jl b/test/testJacobianFoldDeflation.jl
index 7a8a05ee..4fd2e0cc 100644
--- a/test/testJacobianFoldDeflation.jl
+++ b/test/testJacobianFoldDeflation.jl
@@ -151,7 +151,7 @@ res_exp = debugTmpForσ \ rhs
 opt_newton = NewtonPar(tol = 1e-8, verbose = false, eigsolver = EigKrylovKit())
 opts_br0 = ContinuationPar(dsmin = 0.01, dsmax = 0.15, ds= 0.01, p_max = 4.1, max_steps = 250, newton_options = opt_newton, detect_fold = true, detect_bifurcation = 1, nev = 15)
 
-br = continuation(BK.re_make(prob;record_from_solution = (x,p)->norm(x,Inf64)), PALC(), opts_br0, plot = false, verbosity = 0)
+br = continuation(BK.re_make(prob;record_from_solution = (x,p;k...)->norm(x,Inf64)), PALC(), opts_br0, plot = false, verbosity = 0)
 
 opts_br0 = ContinuationPar(dsmin = 0.01, dsmax = 0.15, ds= 0.01, p_max = 4.1, max_steps = 250, newton_options = NewtonPar(tol =1e-8), detect_fold = true, detect_bifurcation = 1, nev = 15)
 
diff --git a/test/testLure.jl b/test/testLure.jl
index 37d13330..2ace5065 100644
--- a/test/testLure.jl
+++ b/test/testLure.jl
@@ -3,7 +3,7 @@ using LinearAlgebra, Test
 using BifurcationKit, Test
 const BK = BifurcationKit
 
-recordFromSolution(x, p) = (u1 = x[1], u2 = x[2])
+recordFromSolution(x, p; k...) = (u1 = x[1], u2 = x[2])
 ####################################################################################################
 function lur!(dz, u, p, t)
     (; α, β) = p
@@ -34,7 +34,7 @@ function plotPO(x, p; k...)
 end
 
 # record function
-function recordPO(x, p)
+function recordPO(x, p; k...)
 	xtt = get_periodic_orbit(p.prob, x, p.p)
 	period = getperiod(p.prob, x, p.p)
 	return (max = maximum(xtt[1,:]), min = minimum(xtt[1,:]), period = period)
@@ -83,7 +83,7 @@ br_po_pd = continuation(br_po, 1, setproperties(br_po.contparams, detect_bifurca
     verbosity = 0, plot = false,
     ampfactor = .2, δp = -0.005,
     usedeflation = true,
-    record_from_solution = (x, p) -> (xtt=reshape(x[1:end-1],3,Mt); return (max = maximum(xtt[1,:]), min = minimum(xtt[1,:]), period = x[end])),
+    record_from_solution = (x, p; k...) -> (xtt=reshape(x[1:end-1],3,Mt); return (max = maximum(xtt[1,:]), min = minimum(xtt[1,:]), period = x[end])),
     normC = norminf
     )
 
diff --git a/test/testNF.jl b/test/testNF.jl
index 61c423c6..daf53252 100644
--- a/test/testNF.jl
+++ b/test/testNF.jl
@@ -84,7 +84,7 @@ br2 = continuation(br, 1, setproperties(opts_br; p_max = 0.2, ds = 0.01, max_ste
 ####################################################################################################
 # Case of the pitchfork
 par_pf = setproperties(prob.params ; x2 = 0.0, x3 = -1.0)
-prob_pf = BK.re_make(prob, params = par_pf, record_from_solution = (x,p)->(x[1], norm(x)))
+prob_pf = BK.re_make(prob, params = par_pf, record_from_solution = (x,p;k...)->(x[1], norm(x)))
 brp = BK.continuation(prob_pf, PALC(tangent=Bordered()), opts_br; normC = norminf)
 
 bpp = BK.get_normal_form(brp, 1; verbose=true)
diff --git a/test/test_bif_detection.jl b/test/test_bif_detection.jl
index 2de18f92..5f4dda74 100644
--- a/test/test_bif_detection.jl
+++ b/test/test_bif_detection.jl
@@ -106,7 +106,7 @@ Jac_m = (x, p; k = 2) -> diagm(0 => p .- x.^k)
 
 opts = ContinuationPar(dsmax = 0.1, dsmin = 1e-5, ds = 0.001, max_steps = 130, p_min = -3., p_max = 0.1, newton_options = NewtonPar(tol = 1e-8, verbose = false, max_iterations = 4), detect_bifurcation=3, n_inversion=4)
 
-prob4 = BK.BifurcationProblem(F, zeros(1), -0.1, (@optic _); J = Jac_m, record_from_solution = (x,p)->x[1])
+prob4 = BK.BifurcationProblem(F, zeros(1), -0.1, (@optic _); J = Jac_m, record_from_solution = (x,p;k...)->x[1])
 br4 = continuation(prob4, alg, opts; verbosity = 0, plot=false)
 testBranch(br4)
 ####################################################################################################
@@ -123,7 +123,7 @@ par = (p1 = -3., p2=-3., k=3)
 
 opts = ContinuationPar(dsmax = 0.1, ds = 0.001, max_steps = 135, p_min = -3., p_max = 4.0, newton_options = NewtonPar(max_iterations = 5), detect_bifurcation = 3, n_inversion = 6, dsmin_bisection = 1e-9, max_bisection_steps = 15, nev = 2)
 
-prob = BK.BifurcationProblem(Ftb, -2ones(2), par, (@optic _.p1); record_from_solution = (x,p)->x[1])
+prob = BK.BifurcationProblem(Ftb, -2ones(2), par, (@optic _.p1); record_from_solution = (x,p;k...)->x[1])
 br = continuation(prob, alg, (@set opts.detect_bifurcation = 3);
     plot = false, verbosity = 0,)
     show(br)