diff --git a/src/periodicorbit/PeriodicOrbitCollocation.jl b/src/periodicorbit/PeriodicOrbitCollocation.jl
index 18dd473b..05df359a 100644
--- a/src/periodicorbit/PeriodicOrbitCollocation.jl
+++ b/src/periodicorbit/PeriodicOrbitCollocation.jl
@@ -677,6 +677,18 @@ function analytical_jacobian_sparse(coll::PeriodicOrbitOCollProblem,
     block_to_sparse(jacBlock)
 end
 
+function jacobian_poocoll_block(coll::PeriodicOrbitOCollProblem,
+                                u::AbstractVector,
+                                pars;
+                                kwargs...) 
+    n, m, Ntst = size(coll)
+    blocks = n * ones(Int64, 1 + m * Ntst + 1); blocks[end] = 1
+    n_blocks = length(blocks)
+    J = BlockArray(zeros(length(u), length(u)), blocks,  blocks)
+    jacobian_poocoll_block!(J, coll, u, pars; kwargs...)
+    return J
+end
+
 @views function jacobian_poocoll_block!(J,
                                 coll::PeriodicOrbitOCollProblem,
                                 u::AbstractVector{𝒯},
@@ -933,7 +945,7 @@ function _newton_pocoll(probPO::PeriodicOrbitOCollProblem,
         _J = analytical_jacobian_sparse(probPO, orbitguess, par)
         jac = (x, p) -> analytical_jacobian!(_J, probPO, x, p)
     else
-        jac = (x, p) -> ForwardDiff.jacobian(z -> probPO(z, p), x)
+        jac = (x, p) -> ForwardDiff.jacobian(z -> residual(probPO, z, p), x)
     end
 
     if options.linsolver isa COPLS
diff --git a/src/periodicorbit/codim2/codim2.jl b/src/periodicorbit/codim2/codim2.jl
index 6e9fd252..500be380 100644
--- a/src/periodicorbit/codim2/codim2.jl
+++ b/src/periodicorbit/codim2/codim2.jl
@@ -291,7 +291,7 @@ function _continuation(gh::Bautin, br::AbstractResult{Tkind, Tprob},
     probsh_fold = BifurcationProblem((x, p) -> residual(pbwrap, x, p), orbitguess, getparams(pbwrap), getlens(pbwrap);
                 J = (x, p) -> jacobian(pbwrap, x, p),
                 Jᵗ = Jᵗ,
-                d2F = (x, p, dx1, dx2) -> d2PO(z -> probPO(z, p), x, dx1, dx2),
+                d2F = (x, p, dx1, dx2) -> d2PO(z -> residual(probPO, z, p), x, dx1, dx2),
                 record_from_solution = _recordsol,
                 plot_solution = _plotsol,
                 )
diff --git a/test/cop.jl b/test/cop.jl
index dfbcf220..50aaa7b7 100644
--- a/test/cop.jl
+++ b/test/cop.jl
@@ -20,7 +20,7 @@ prob_col = PeriodicOrbitOCollProblem(Ntst, m;
                                     xπ = rand(N*( 1 + m * Ntst)))
 _ci = generate_solution(prob_col, t->cos(t) .* ones(N), 2pi);
 #####################################################
-Jcofd = ForwardDiff.jacobian(z->prob_col(z, par_sl), _ci);
+Jcofd = ForwardDiff.jacobian(z->BK.residual(prob_col,z, par_sl), _ci);
 Jco = BK.analytical_jacobian(prob_col, _ci, par_sl);
 
 _rhs = rand(size(Jco, 1))
diff --git a/test/stuartLandauCollocation.jl b/test/stuartLandauCollocation.jl
index 76ae099d..63166ea8 100644
--- a/test/stuartLandauCollocation.jl
+++ b/test/stuartLandauCollocation.jl
@@ -52,7 +52,7 @@ _ci = BK.generate_solution(prob_col, _orbit, 1.)
 BK.get_periodic_orbit(prob_col, _ci, par_sl)
 BK.getmaximum(prob_col, _ci, par_sl)
 @test BK.∂(sin, Val(2))(0.) == 0
-prob_col(_ci, par_sl) #|> scatter
+BK.residual(prob_col, _ci, par_sl) #|> scatter
 BK.get_time_slices(prob_col, _ci)
 
 # interpolate solution
@@ -156,13 +156,13 @@ prob_col.ϕ[2] = 1 #phase condition
 
 _orbit(t) = [cos(t), sin(t)] * sqrt(par_sl.r/par_sl.c3)
 _ci = BK.generate_solution(prob_col, _orbit, 2pi)
-prob_col(_ci, par_sl)
-@test prob_col(_ci, par_sl)[1:end-1] |> norminf < 1e-7
+BK.residual(prob_col, _ci, par_sl)
+@test BK.residual(prob_col, _ci, par_sl)[1:end-1] |> norminf < 1e-7
 
 prob_coll_ip = @set prob_col.prob_vf = probsl_ip
 
-@time prob_col(_ci, par_sl);
-@time prob_coll_ip(_ci, par_sl);
+@time BK.residual(prob_col, _ci, par_sl);
+@time BK.residual(prob_coll_ip, _ci, par_sl);
 
 # test precision of generated solution
 _sol = BK.get_periodic_orbit(prob_col, _ci, nothing)
@@ -222,8 +222,8 @@ nullvf(x,p) = zero(x)
 prob0 = BifurcationProblem(nullvf, zeros(N), par_hopf, (@optic _.r))
 prob_col = BK.PeriodicOrbitOCollProblem(Ntst, m; prob_vf = prob0, N = N, ϕ = rand(N*( 1 + m * Ntst)), xπ = rand(N*( 1 + m * Ntst)))
 _ci = BK.generate_solution(prob_col, t->cos(t) .* ones(N), 2pi);
-prob_col(_ci, par_sl);
-Jcofd = ForwardDiff.jacobian(z->prob_col(z, par_sl), _ci);
+BK.residual(prob_col,_ci, par_sl);
+Jcofd = ForwardDiff.jacobian(z -> BK.residual(prob_col, z, par_sl), _ci);
 D = @time BK.analytical_jacobian(prob_col, _ci, par_sl); #0.000121 seconds (341 allocations: 156.516 KiB)
 @test norminf(Jcofd - D) < 1e-15
 
@@ -236,7 +236,7 @@ idvf(x,p) = _al*x
 prob_ana = BifurcationProblem(idvf, zeros(N), par_hopf, (@optic _.r) ; J = (x,p) -> _al)
 prob_col = BK.PeriodicOrbitOCollProblem(Ntst, m; prob_vf = prob_ana, N = N, ϕ = rand(N*( 1 + m * Ntst)), xπ = rand(N*( 1 + m * Ntst)))
 _ci = BK.generate_solution(prob_col, t->cos(t) .* ones(N), 2pi);
-Jcofd = ForwardDiff.jacobian(z -> prob_col(z, par_sl), _ci);
+Jcofd = ForwardDiff.jacobian(z -> BK.residual(prob_col, z, par_sl), _ci);
 Jco = BK.analytical_jacobian(prob_col, _ci, par_sl); # 0.004388 seconds (573 allocations: 60.124 MiB)
 @test norminf(Jcofd - Jco) < 1e-15
 
@@ -247,7 +247,7 @@ Ntst = 3
 m = 2
 prob_col = BK.PeriodicOrbitOCollProblem(Ntst, m; prob_vf = probsl, N = N, ϕ = rand(N*( 1 + m * Ntst)), xπ = rand(N*( 1 + m * Ntst)))
 _ci = BK.generate_solution(prob_col, t->cos(t) .* ones(N), 2pi);
-Jcofd = ForwardDiff.jacobian(z->prob_col(z, par_sl), _ci);
+Jcofd = ForwardDiff.jacobian(z->BK.residual(prob_col, z, par_sl), _ci);
 Jco = @time BK.analytical_jacobian(prob_col, _ci, par_sl);
 Jco_bk = @time BK.jacobian_poocoll_block(prob_col, _ci, par_sl);
 @test norminf(Jcofd - Jco) < 1e-14
diff --git a/test/testHopfMA.jl b/test/testHopfMA.jl
index 43f96f9a..510c67f9 100644
--- a/test/testHopfMA.jl
+++ b/test/testHopfMA.jl
@@ -204,7 +204,7 @@ prob = BifurcationKit.BifurcationProblem(Fbru!, sol0, par_bru, (@optic _.l);
 
 poTrap = PeriodicOrbitTrapProblem(prob, real.(vec_hopf), hopfpt.u, M, 2n    )
 
-jac_PO_fd = BK.finite_differences(x -> poTrap(x, (@set par_bru.l = l_hopf + 0.01)), orbitguess_f)
+jac_PO_fd = BK.finite_differences(x -> BK.residual(poTrap, x, (@set par_bru.l = l_hopf + 0.01)), orbitguess_f)
 jac_PO_sp = poTrap(Val(:JacFullSparse), orbitguess_f, (@set par_bru.l = l_hopf + 0.01))
 
 # test of the Jacobian for PeriodicOrbit via Finite differences VS the FD associated jacobian
@@ -218,7 +218,7 @@ BK.get_time_diff(poTrap, orbitguess_f)
 # BK.Jc(poTrap, orbitguess_f, par_bru, orbitguess_f)
 
 # newton to find Periodic orbit
-_prob = BK.BifurcationProblem((x, p) -> poTrap(x, p), copy(orbitguess_f), (@set par_bru.l = l_hopf + 0.01); J = (x, p) ->  poTrap(Val(:JacFullSparse),x,p))
+_prob = BK.BifurcationProblem((x, p) -> BK.residual(poTrap, x, p), copy(orbitguess_f), (@set par_bru.l = l_hopf + 0.01); J = (x, p) ->  poTrap(Val(:JacFullSparse),x,p))
 opt_po = NewtonPar(tol = 1e-8, max_iterations = 150)
 outpo_f = @time BK.solve(_prob, Newton(), opt_po)
 @test BK.converged(outpo_f)
diff --git a/test/test_potrap.jl b/test/test_potrap.jl
index 48fd62e1..27e496b2 100644
--- a/test/test_potrap.jl
+++ b/test/test_potrap.jl
@@ -31,9 +31,9 @@ pbi = PeriodicOrbitTrapProblem(
 @test BK.isinplace(pb) == false
 # BK.POTrapFunctional(pb, res, orbitguess_f)
 # BK.POTrapFunctional(pbi, res, orbitguess_f)
-res = pb(orbitguess_f, par)
-resg = pbg(orbitguess_f, par)
-resi = pbi(orbitguess_f, par)
+res = BK.residual(pb, orbitguess_f, par)
+resg = BK.residual(pbg, orbitguess_f, par)
+resi = BK.residual(pbi, orbitguess_f, par)
 @test res == resi
 @test res == resg
 
@@ -43,7 +43,7 @@ resi = pbi(orbitguess_f, par, orbitguess_f)
 @test res == resi
 @test res == resg
 
-BK.potrap_functional!(pbi, resi, orbitguess_f, par)
+BK.residual!(pbi, resi, orbitguess_f, par)
 BK.potrap_functional_jac!(pbi, resi, orbitguess_f, par, orbitguess_f)
 @test res == resi
 
@@ -147,7 +147,7 @@ function _dfunctional(poPb, u0, p, du)
 end
 
 
-res = pb(orbitguess_f, par)
+res = BK.residual(pb, orbitguess_f, par)
 _res = _functional(pb, orbitguess_f, par)
 @test res ≈ _res
 
@@ -158,7 +158,7 @@ _res = _dfunctional(pb, orbitguess_f, par, _du)
 
 # with mass matrix
 pbmass = @set pb.massmatrix = spdiagm( 0 => rand(pb.N))
-res = pbmass(orbitguess_f, par)
+res = BK.residual(pbmass, orbitguess_f, par)
 _res = _functional(pbmass, orbitguess_f, par)
 @test res ≈ _res
 
@@ -177,7 +177,7 @@ pbsp = PeriodicOrbitTrapProblem(
             10)
 orbitguess_f = rand(2n*10+1)
 dorbit = rand(2n*10+1)
-Jfd = sparse( ForwardDiff.jacobian(x -> pbsp(x, par), orbitguess_f) )
+Jfd = sparse( ForwardDiff.jacobian(x -> BK.residual(pbsp,x, par), orbitguess_f) )
 Jan = pbsp(Val(:JacFullSparse), orbitguess_f, par)
 @test norm(Jfd - Jan, Inf) < 1e-6
 
@@ -199,7 +199,7 @@ Jan = pbsp(Val(:JacCyclicSparse), orbitguess_f, par)
 
 # we update the section
 BK.updatesection!(pbsp, rand(2n*10+1), par)
-Jfd2 = sparse( ForwardDiff.jacobian(x -> pbsp(x, par), orbitguess_f) )
+Jfd2 = sparse( ForwardDiff.jacobian(x -> BK.residual(pbsp, x, par), orbitguess_f) )
 Jan2 = pbsp(Val(:JacFullSparse), orbitguess_f, par)
 @test Jan2 != Jan
 @test norm(Jfd2 - Jan2, Inf) < 1e-6
@@ -211,7 +211,7 @@ Jan = pbsp(Val(:JacCyclicSparse), orbitguess_f, par)
 ##########################
 #### idem with mass matrix
 pbsp_mass = @set pbsp.massmatrix = massmatrix = spdiagm( 0 => rand(pbsp.N))
-Jfd = sparse( ForwardDiff.jacobian(x -> pbsp_mass(x, par), orbitguess_f) )
+Jfd = sparse( ForwardDiff.jacobian(x -> BK.residual(pbsp_mass, x, par), orbitguess_f) )
 Jan = pbsp_mass(Val(:JacFullSparse), orbitguess_f, par)
 @test norm(Jfd - Jan, Inf) < 1e-6
 
@@ -229,11 +229,11 @@ pbsp_mass(Val(:JacFullSparseInplace), Jan2, orbitguess_f, par, _indx)
 @test norm(Jan2 - Jan, Inf) < 1e-6
 
 Jan = pbsp_mass(Val(:JacCyclicSparse), orbitguess_f, par)
-@test norm(Jfd[1:size(Jan,1),1:size(Jan,1)] - Jan, Inf) < 1e-6
+@test norm(Jfd[1:size(Jan,1), 1:size(Jan,1)] - Jan, Inf) < 1e-6
 
 # we update the section
 BK.updatesection!(pbsp_mass, rand(2n*10+1), par)
-Jfd2 = sparse( ForwardDiff.jacobian(x -> pbsp_mass(x, par), orbitguess_f) )
+Jfd2 = sparse( ForwardDiff.jacobian(x -> BK.residual(pbsp_mass, x, par), orbitguess_f) )
 Jan2 = pbsp_mass(Val(:JacFullSparse), orbitguess_f, par)
 @test Jan2 != Jan
 @test norm(Jfd2 - Jan2, Inf) < 1e-6
@@ -315,5 +315,5 @@ BK.get_periodic_orbit(pbspti, orbitguess_f, par)
 
 @test pbspti.xπ ≈ pbsp.xπ
 @test pbspti.ϕ ≈ pbsp.ϕ
-pbspti(orbitguess_f, par)
-@test pbsp(orbitguess_f, par) ≈ pbspti(orbitguess_f, par)
+BK.residual(pbspti, orbitguess_f, par)
+@test BK.residual(pbsp, orbitguess_f, par) ≈ BK.residual(pbspti, orbitguess_f, par)
\ No newline at end of file