diff --git a/Project.toml b/Project.toml
index 33270f02..1259a0db 100644
--- a/Project.toml
+++ b/Project.toml
@@ -11,9 +11,9 @@ CommonSolve = "38540f10-b2f7-11e9-35d8-d573e4eb0ff2"
 DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae"
 
 [compat]
-CTBase = "0.12"
-CTDirect = "0.11"
-CTFlows = "0.5"
+CTBase = "0.13"
+CTDirect = "0.12"
+CTFlows = "0.6"
 CommonSolve = "0.2"
 DocStringExtensions = "0.9"
 julia = "1.10"
diff --git a/docs/Project.toml b/docs/Project.toml
index f6e374ac..ae87e81c 100644
--- a/docs/Project.toml
+++ b/docs/Project.toml
@@ -22,5 +22,22 @@ Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
 Roots = "f2b01f46-fcfa-551c-844a-d8ac1e96c665"
 
 [compat]
+BenchmarkTools = "1.5"
+CTBase = "0.13"
+CTDirect = "0.12"
+CTFlows = "0.6"
+CommonSolve = "0.2"
+DifferentiationInterface = "0.5"
+Documenter = "1.6"
+DocumenterMermaid = "0.1"
+ForwardDiff = "0.10"
+Interpolations = "0.15"
+JLD2 = "0.4"
+JSON3 = "1.14"
 MINPACK = "1.3"
+MadNLP = "0.8"
+NLPModelsIpopt = "0.10"
+Percival = "0.7"
+Plots = "1.40"
+Roots = "2.1"
 julia = "1.10"
diff --git a/docs/src/tutorial-basic-example.md b/docs/src/tutorial-basic-example.md
index 67ab5bb2..20655e76 100644
--- a/docs/src/tutorial-basic-example.md
+++ b/docs/src/tutorial-basic-example.md
@@ -113,7 +113,7 @@ println("Objective from loaded solution: ", sol_reloaded.objective)
 
 # JSON export / read
 export_ocp_solution(sol, filename_prefix="my_solution")
-sol_json = import_ocp_solution("my_solution")
+sol_json = import_ocp_solution(ocp, filename_prefix="my_solution")
 println("Objective from JSON discrete solution: ", sol_json.objective)
 ```
 
diff --git a/docs/src/tutorial-continuation.md b/docs/src/tutorial-continuation.md
index 6b28f8af..7d330eed 100644
--- a/docs/src/tutorial-continuation.md
+++ b/docs/src/tutorial-continuation.md
@@ -103,7 +103,7 @@ objective!(ocp, :mayer, (x0, xf, v) -> xf[1], :max)
 dynamics!(ocp, (x, u, v) -> F0(x) + u*F1(x) )
 
 sol0 = solve(ocp; display=false)
-@printf("Objective for reference solution %.6f\n", sol0.objective)
+@printf("Objective for reference solution %.6f\n", objective(sol0))
 ```
 
 Then we perform the continuation on the maximal thrust.
@@ -115,9 +115,9 @@ obj_list  = []
 for Tmax_local=3.5:-0.5:1
     global Tmax = Tmax_local  
     global sol = solve(ocp; display=false, init=sol)
-    @printf("Tmax %.2f objective %.6f iterations %d\n", Tmax, sol.objective, sol.iterations)
+    @printf("Tmax %.2f objective %.6f iterations %d\n", Tmax, objective(sol), iterations(sol))
     push!(Tmax_list, Tmax)
-    push!(obj_list, sol.objective)
+    push!(obj_list, objective(sol))
 end 
 ```
 
diff --git a/docs/src/tutorial-flow.md b/docs/src/tutorial-flow.md
index 0cb790ea..71f1f6d0 100644
--- a/docs/src/tutorial-flow.md
+++ b/docs/src/tutorial-flow.md
@@ -139,7 +139,7 @@ plot(sol)
 You can notice from the graph of `v` that the integrator has made very few steps:
 
 ```@example main
-sol.times
+time_grid(sol)
 ```
 
 To have a better visualisation (the accuracy won't change), you can provide a fine grid.
diff --git a/docs/src/tutorial-goddard.md b/docs/src/tutorial-goddard.md
index a8c20249..8445a407 100644
--- a/docs/src/tutorial-goddard.md
+++ b/docs/src/tutorial-goddard.md
@@ -124,10 +124,10 @@ arc is with zero control. Note that the switching function vanishes along the si
 boundary arcs.
 
 ```@example main
-t = direct_sol.times
-x = direct_sol.state
-u = direct_sol.control
-p = direct_sol.costate
+t = time_grid(direct_sol)
+x = state(direct_sol)
+u = control(direct_sol)
+p = costate(direct_sol)
 
 H1 = Lift(F1)           # H1(x, p) = p' * F1(x)
 φ(t) = H1(x(t), p(t))   # switching function
diff --git a/docs/src/tutorial-initial-guess.md b/docs/src/tutorial-initial-guess.md
index 82881c8e..7370c17b 100644
--- a/docs/src/tutorial-initial-guess.md
+++ b/docs/src/tutorial-initial-guess.md
@@ -43,7 +43,7 @@ This will default to initialize all variables to 0.1.
 ```@example main
 # solve the optimal control problem without initial guess
 sol = solve(ocp; display=false)
-println("Number of iterations: ", sol.iterations)
+println("Number of iterations: ", iterations(sol))
 nothing # hide
 ```
 
@@ -57,10 +57,10 @@ Note that the following formulations are equivalent to not giving an initial gue
 
 ```@example main
 sol = solve(ocp; init=nothing, display=false)
-println("Number of iterations: ", sol.iterations)
+println("Number of iterations: ", iterations(sol))
 
 sol = solve(ocp; init=(), display=false)
-println("Number of iterations: ", sol.iterations)
+println("Number of iterations: ", iterations(sol))
 nothing # hide
 ```
 
@@ -74,7 +74,7 @@ We first illustrate the constant initial guess, using vectors or scalars accordi
 ```@example main
 # solve the optimal control problem with initial guess with constant values
 sol = solve(ocp; init=(state=[-0.2, 0.1], control=-0.2, variable=0.05), display=false)
-println("Number of iterations: ", sol.iterations)
+println("Number of iterations: ", iterations(sol))
 nothing # hide
 ```
 
@@ -82,15 +82,15 @@ Partial initializations are also valid, as shown below. Note the ending comma wh
 ```@example main
 # initialisation only on the state
 sol = solve(ocp; init=(state=[-0.2, 0.1],), display=false)
-println("Number of iterations: ", sol.iterations)
+println("Number of iterations: ", iterations(sol))
 
 # initialisation only on the control
 sol = solve(ocp; init=(control=-0.2,), display=false)
-println("Number of iterations: ", sol.iterations)
+println("Number of iterations: ", iterations(sol))
 
 # initialisation only on the state and the variable
 sol = solve(ocp; init=(state=[-0.2, 0.1], variable=0.05), display=false)
-println("Number of iterations: ", sol.iterations)
+println("Number of iterations: ", iterations(sol))
 nothing # hide
 ```
 
@@ -103,7 +103,7 @@ x(t) = [ -0.2t, 0.1t ]
 u(t) = -0.2t
 
 sol = solve(ocp; init=(state=x, control=u, variable=0.05), display=false)
-println("Number of iterations: ", sol.iterations)
+println("Number of iterations: ", iterations(sol))
 nothing # hide
 ```
 
@@ -120,7 +120,7 @@ x_vec = [[0, 0], [-0.1, 0.3], [-0.15,0.4], [-0.3, 0.5]]
 u_vec = [0, -0.8,  -0.3, 0]
 
 sol = solve(ocp; init=(time=t_vec, state=x_vec, control=u_vec, variable=0.05), display=false)
-println("Number of iterations: ", sol.iterations)
+println("Number of iterations: ", iterations(sol))
 nothing # hide
 ```
 
@@ -133,11 +133,11 @@ The constant, functional and vector initializations can be mixed, for instance a
 ```@example main
 # we can mix constant values with functions of time
 sol = solve(ocp; init=(state=[-0.2, 0.1], control=u, variable=0.05), display=false)
-println("Number of iterations: ", sol.iterations)
+println("Number of iterations: ", iterations(sol))
 
 # wa can mix every possibility
 sol = solve(ocp; init=(time=t_vec, state=x_vec, control=u, variable=0.05), display=false)
-println("Number of iterations: ", sol.iterations)
+println("Number of iterations: ", iterations(sol))
 nothing # hide
 ```
 
@@ -153,24 +153,24 @@ sol_init = solve(ocp; display=false)
 
 # solve the problem using solution as initial guess
 sol = solve(ocp; init=sol_init, display=false)
-println("Number of iterations: ", sol.iterations)
+println("Number of iterations: ", iterations(sol))
 nothing # hide
 ```
 
 Note that you can also manually pick and choose which data to reuse from a solution, by recovering the 
-functions ```sol.state```, ```sol.control``` and the values ```sol.variable```.
+functions ```state(sol)```, ```control(sol)``` and the values ```variable(sol)```.
 For instance the following formulation is equivalent to the ```init=sol``` one.
 
 ```@example main
 # use a previous solution to initialise picking data
 sol = solve(ocp; 
     init = (
-        state    = sol.state, 
-        control  = sol.control, 
-        variable = sol.variable
+        state    = state(sol), 
+        control  = control(sol), 
+        variable = variable(sol)
     ), 
     display=false)
-println("Number of iterations: ", sol.iterations)
+println("Number of iterations: ", iterations(sol))
 nothing # hide
 ``` 
 
diff --git a/docs/src/tutorial-iss.md b/docs/src/tutorial-iss.md
index 164bc958..2b4e8545 100644
--- a/docs/src/tutorial-iss.md
+++ b/docs/src/tutorial-iss.md
@@ -284,8 +284,8 @@ function pretty_plot(S, p0; Np0=20, kwargs...)
     p0s = range(p0_min, p0_max, length=Np0)
     for i ∈ eachindex(p0s)
         sol = exp(p0s[i])
-        x = [sol.state(t)   for t ∈ sol.times]
-        p = [sol.costate(t) for t ∈ sol.times]
+        x = [state(sol)(t)   for t ∈ time_grid(sol)]
+        p = [costate(sol)(t) for t ∈ time_grid(sol)]
         label = i==1 ? "extremals" : false
         plot!(plt_flow, x, p, color=:blue, label=label)
     end
@@ -296,8 +296,8 @@ function pretty_plot(S, p0; Np0=20, kwargs...)
     ps  = zeros(length(p0s), length(times))
     for i ∈ eachindex(p0s)
         sol = exp(p0s[i], saveat=times)
-        xs[i, :] .= sol.state.(times)
-        ps[i, :] .= sol.costate.(times)
+        xs[i, :] .= state(sol).(times)
+        ps[i, :] .= costate(sol).(times)
     end
     for j ∈ eachindex(times)
         label = j==1 ? "flow at times" : false
@@ -310,8 +310,8 @@ function pretty_plot(S, p0; Np0=20, kwargs...)
     
     # solution
     sol = exp(p0_sol)
-    x = [sol.state(t)   for t ∈ sol.times]
-    p = [sol.costate(t) for t ∈ sol.times]
+    x = [state(sol)(t)   for t ∈ time_grid(sol)]
+    p = [costate(sol)(t) for t ∈ time_grid(sol)]
     plot!(plt_flow, x, p, color=:red, linewidth=2, label="extremal solution")
     plot!(plt_flow, [x[end]], [p[end]], seriestype=:scatter, color=:green, label=false)
 
@@ -347,8 +347,8 @@ function pretty_plot(S, p0; Np0=20, kwargs...) # hide
     p0s = range(p0_min, p0_max, length=Np0) # hide
     for i ∈ eachindex(p0s) # hide
         sol = exp(p0s[i]) # hide
-        x = [sol.state(t)   for t ∈ sol.times] # hide
-        p = [sol.costate(t) for t ∈ sol.times] # hide
+        x = [state(sol)(t)   for t ∈ time_grid(sol)] # hide
+        p = [costate(sol)(t) for t ∈ time_grid(sol)] # hide
         label = i==1 ? "extremals" : false # hide
         plot!(plt_flow, x, p, color=:blue, label=label) # hide
     end # hide
@@ -359,8 +359,8 @@ function pretty_plot(S, p0; Np0=20, kwargs...) # hide
     ps  = zeros(length(p0s), length(times)) # hide
     for i ∈ eachindex(p0s) # hide
         sol = exp(p0s[i], saveat=times) # hide
-        xs[i, :] .= sol.state.(times) # hide
-        ps[i, :] .= sol.costate.(times) # hide
+        xs[i, :] .= state(sol).(times) # hide
+        ps[i, :] .= costate(sol).(times) # hide
     end # hide
     for j ∈ eachindex(times) # hide
         label = j==1 ? "flow at times" : false # hide
@@ -373,8 +373,8 @@ function pretty_plot(S, p0; Np0=20, kwargs...) # hide
      # hide
     # solution # hide
     sol = exp(p0_sol) # hide
-    x = [sol.state(t)   for t ∈ sol.times] # hide
-    p = [sol.costate(t) for t ∈ sol.times] # hide
+    x = [state(sol)(t)   for t ∈ time_grid(sol)] # hide
+    p = [costate(sol)(t) for t ∈ time_grid(sol)] # hide
     plot!(plt_flow, x, p, color=:red, linewidth=2, label="extremal solution") # hide
     plot!(plt_flow, [x[end]], [p[end]], seriestype=:scatter, color=:green, label=false) # hide
  # hide
diff --git a/docs/src/tutorial-plot.md b/docs/src/tutorial-plot.md
index 4b2ed6be..3dac0acc 100644
--- a/docs/src/tutorial-plot.md
+++ b/docs/src/tutorial-plot.md
@@ -79,7 +79,7 @@ optimal control problem.
 t0 = 0
 tf = 1
 x0 = [ -1, 0 ]
-p0 = sol.costate(t0)
+p0 = costate(sol)(t0)
 f  = Flow(ocp, (x, p) -> p[2])
 sol_flow = f( (t0, tf), x0, p0 )
 plot(sol_flow)
@@ -152,10 +152,10 @@ You can of course create your own plots by getting the `state`, `costate` and `c
 
 ```@example main
 using LinearAlgebra
-t = sol.times
-x = sol.state
-p = sol.costate
-u = sol.control
+t = time_grid(sol)
+x = state(sol)
+p = costate(sol)
+u = control(sol)
 plot(t, norm∘u; label="‖u‖") 
 ```
 
@@ -163,7 +163,7 @@ plot(t, norm∘u; label="‖u‖")
 
     - The `norm` function is from `LinearAlgebra.jl`. 
     - The `∘` operator is the composition operator. Hence, `norm∘u` is the function `t -> norm(u(t))`. 
-    - The `sol.state`, `sol.costate` and `sol.control` are functions that return the state, costate and control trajectories at a given time.
+    - The `state(sol)`, `costate(sol)` and `control(sol)` are functions that return the state, costate and control trajectories at a given time.
 
 
 ## Normalized time
diff --git a/docs/src/tutorial-solve.md b/docs/src/tutorial-solve.md
index cb2bddd8..298df3ff 100644
--- a/docs/src/tutorial-solve.md
+++ b/docs/src/tutorial-solve.md
@@ -153,14 +153,14 @@ There are `init`, `grid_size` and `time_grid`.
 
 ```@example main
 sol = solve(ocp; grid_size=10, display=false)
-sol.times
+time_grid(sol)
 ```
 
 Or with MadNLP.jl:
 
 ```@example main
 sol = solve(ocp, :madnlp; grid_size=10, display=false)
-sol.times
+time_grid(sol)
 ```
 
 ### NLPModelsIpopt
diff --git a/src/OptimalControl.jl b/src/OptimalControl.jl
index d49c55d8..89ee0fdb 100644
--- a/src/OptimalControl.jl
+++ b/src/OptimalControl.jl
@@ -53,9 +53,10 @@ export OptimalControlModel, OptimalControlSolution
 export Autonomous, NonAutonomous
 export NonFixed, Fixed
 export Model, __OCPModel
-export variable!,
-    time!, constraint!, dynamics!, objective!, state!, control!, remove_constraint!, constraint
-#export is_time_independent, is_time_dependent, is_min, is_max, is_variable_dependent, is_variable_independent
+export variable!, time!, constraint!, dynamics!, objective!, state!, control!, remove_constraint!
+export constraint
+export time_grid, control, state, variable, costate, objective
+export iterations, stopping, message, infos
 export Lie, @Lie, Poisson, Lift, ⋅, ∂ₜ
 export @def
 export ct_repl, ct_repl_update_model
diff --git a/test/Project.toml b/test/Project.toml
index 961bd124..20a42acb 100644
--- a/test/Project.toml
+++ b/test/Project.toml
@@ -8,4 +8,9 @@ SciMLBase = "0bca4576-84f4-4d90-8ffe-ffa030f20462"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [compat]
+CTDirect = "0.12"
+NLPModelsIpopt = "0.10"
+NonlinearSolve = "3.14"
+OrdinaryDiffEq = "6.88"
+SciMLBase = "2.48"
 julia = "1.10"
diff --git a/test/test_basic.jl b/test/test_basic.jl
index 2d8430dd..efff715f 100644
--- a/test/test_basic.jl
+++ b/test/test_basic.jl
@@ -10,5 +10,5 @@ function test_basic()
     end
 
     sol = solve(ocp; display = false)
-    @test sol.objective ≈ 6 atol = 1e-2
+    @test objective(sol) ≈ 6 atol = 1e-2
 end
diff --git a/test/test_continuation.jl b/test/test_continuation.jl
index 0a1b5f45..ac7a6288 100644
--- a/test/test_continuation.jl
+++ b/test/test_continuation.jl
@@ -30,7 +30,7 @@ function test_continuation()
                 ocp = ocp_T(T)
                 sol = solve(ocp, print_level = 0, init = init)
                 init = sol
-                push!(obj_list, sol.objective)
+                push!(obj_list, objective(sol))
             end
             @test obj_list ≈ [12, 1.5, 0.44, 0.19, 0.096] rtol = 1e-2
         end
@@ -69,7 +69,7 @@ function test_continuation()
                 ocp = myocp(ρ)
                 sol = solve(ocp, print_level = 0, init = init)
                 init = sol
-                push!(obj_list, sol.objective)
+                push!(obj_list, objective(sol))
             end
             @test obj_list ≈ [-0.034, -1.7, -6.2, -35, -148] rtol = 1e-2
         end
@@ -86,7 +86,7 @@ function test_continuation()
             for Tmax = 3.5:-0.5:1
                 sol = solve(goddard(Tmax = Tmax).ocp, print_level = 0, init = sol)
                 push!(Tmax_list, Tmax)
-                push!(obj_list, sol.objective)
+                push!(obj_list, objective(sol))
             end
             @test obj_list ≈ [1.0125, 1.0124, 1.0120, 1.0112, 1.0092, 1.0036] rtol = 1e-2
 
diff --git a/test/test_goddard_direct.jl b/test/test_goddard_direct.jl
index 130a4048..056d9a9c 100644
--- a/test/test_goddard_direct.jl
+++ b/test/test_goddard_direct.jl
@@ -1,10 +1,10 @@
 function test_goddard_direct()
 
     # goddard with state constraint - maximize altitude
-    ocp, objective = Goddard()
+    ocp, obj = Goddard()
 
     # initial guess (constant state and control functions)
     init = (state = [1.01, 0.05, 0.8], control = 0.1, variable = 0.2)
     sol = solve(ocp; grid_size = 10, display = false, init = init)
-    @test sol.objective ≈ objective atol = 5e-3
+    @test objective(sol) ≈ obj atol = 5e-3
 end
diff --git a/test/test_grid.jl b/test/test_grid.jl
index 0041754c..e58f6aad 100644
--- a/test/test_grid.jl
+++ b/test/test_grid.jl
@@ -13,7 +13,7 @@ function test_grid()
 
     time_grid = [0, 0.1, 0.3, 0.6, 0.98, 0.99, 1]
     sol = solve(ocp; time_grid = time_grid, display = false)
-    @test sol.objective ≈ 0.309 rtol = 1e-2
+    @test objective(sol) ≈ 0.309 rtol = 1e-2
 
     # 1. simple integrator min energy (dual control for test)
     ocp = Model()
@@ -31,14 +31,14 @@ function test_grid()
     @testset verbose = true showtiming = true ":explicit_grid" begin
         time_grid = LinRange(0, 1, CTDirect.__grid_size() + 1)
         sol = solve(ocp, time_grid = time_grid, print_level = 0)
-        @test sol.objective == sol0.objective
-        @test sol.iterations == sol0.iterations
+        @test objective(sol) == objective(sol0)
+        @test iterations(sol) == iterations(sol0)
     end
 
     @testset verbose = true showtiming = true ":non_uniform_grid" begin
         time_grid = [0, 0.1, 0.3, 0.6, 0.98, 0.99, 1]
         sol = solve(ocp, time_grid = time_grid, print_level = 0)
-        @test sol.objective ≈ 0.309 rtol = 1e-2
+        @test objective(sol) ≈ 0.309 rtol = 1e-2
     end
 
     # 2. integrator free times
@@ -58,14 +58,14 @@ function test_grid()
 
     @testset verbose = true showtiming = true ":explicit_grid" begin
         sol = solve(ocp, time_grid = LinRange(0, 1, CTDirect.__grid_size() + 1), print_level = 0)
-        @test sol.objective == sol0.objective
-        @test sol.iterations == sol0.iterations
+        @test objective(sol) == objective(sol0)
+        @test iterations(sol) == iterations(sol0)
     end
 
     @testset verbose = true showtiming = true ":non_uniform_grid" begin
         sol = solve(ocp, time_grid = [0, 0.1, 0.3, 0.5, 0.6, 0.8, 0.95, 1], print_level = 0)
         #plot(sol, show=true)
-        @test sol.objective ≈ 7.96 rtol = 1e-2
+        @test objective(sol) ≈ 7.96 rtol = 1e-2
     end
 
     # 3. parametric ocp (T=2) with explicit / non-uniform grid
@@ -86,12 +86,12 @@ function test_grid()
 
     @testset verbose = true showtiming = true ":explicit_grid" begin
         sol = solve(ocpT2, time_grid = LinRange(0, 1, CTDirect.__grid_size() + 1), print_level = 0)
-        @test sol.objective == sol0.objective
-        @test sol.iterations == sol0.iterations
+        @test objective(sol) == objective(sol0)
+        @test iterations(sol) == iterations(sol0)
     end
 
     @testset verbose = true showtiming = true ":non_uniform_grid" begin
         sol = solve(ocpT2, time_grid = [0, 0.3, 1, 1.9, 2], print_level = 0)
-        @test sol.objective ≈ 2.43 rtol = 1e-2
+        @test objective(sol) ≈ 2.43 rtol = 1e-2
     end
 end
diff --git a/test/test_initial_guess.jl b/test/test_initial_guess.jl
index 0410e790..b3cf45f6 100644
--- a/test/test_initial_guess.jl
+++ b/test/test_initial_guess.jl
@@ -6,13 +6,13 @@ function test_initial_guess()
 
     # test functions
     function check_xf(sol, xf)
-        return xf == sol.state(sol.times[end])
+        return xf == state(sol)(time_grid(sol)[end])
     end
     function check_uf(sol, uf)
-        return uf == sol.control(sol.times[end])
+        return uf == control(sol)(time_grid(sol)[end])
     end
     function check_v(sol, v)
-        return v == sol.variable
+        return v == variable(sol)
     end
 
     # reference solution
@@ -114,11 +114,11 @@ function test_initial_guess()
     # 1. functional initial guess
     @testset verbose = true showtiming = true ":functional_x" begin
         sol = solve(ocp, print_level = 0, init = (state = x_func,), max_iter = maxiter)
-        @test(check_xf(sol, x_func(sol.times[end])))
+        @test(check_xf(sol, x_func(time_grid(sol)[end])))
     end
     @testset verbose = true showtiming = true ":functional_u" begin
         sol = solve(ocp, print_level = 0, init = (control = u_func,), max_iter = maxiter)
-        @test(check_uf(sol, u_func(sol.times[end])))
+        @test(check_uf(sol, u_func(time_grid(sol)[end])))
     end
     @testset verbose = true showtiming = true ":functional_xu" begin
         sol = solve(
@@ -127,7 +127,7 @@ function test_initial_guess()
             init = (state = x_func, control = u_func),
             max_iter = maxiter,
         )
-        @test(check_xf(sol, x_func(sol.times[end])) && check_uf(sol, u_func(sol.times[end])))
+        @test(check_xf(sol, x_func(time_grid(sol)[end])) && check_uf(sol, u_func(time_grid(sol)[end])))
     end
 
     # 1.d interpolated initial guess
@@ -171,7 +171,7 @@ function test_initial_guess()
         )
         @test(
             check_xf(sol, x_vec[end]) &&
-            check_uf(sol, u_func(sol.times[end])) &&
+            check_uf(sol, u_func(time_grid(sol)[end])) &&
             check_v(sol, v_const)
         )
     end
@@ -180,9 +180,9 @@ function test_initial_guess()
     @testset verbose = true showtiming = true ":warm_start" begin
         sol = solve(ocp, print_level = 0, init = sol0, max_iter = maxiter)
         @test(
-            check_xf(sol, sol.state(sol.times[end])) &&
-            check_uf(sol, sol.control(sol.times[end])) &&
-            check_v(sol, sol.variable)
+            check_xf(sol, state(sol)(time_grid(sol)[end])) &&
+            check_uf(sol, control(sol)(time_grid(sol)[end])) &&
+            check_v(sol, variable(sol))
         )
     end
 end
diff --git a/test/test_objective.jl b/test/test_objective.jl
index cef14c74..79cce76a 100644
--- a/test/test_objective.jl
+++ b/test/test_objective.jl
@@ -16,7 +16,7 @@ function test_objective()
 
     @testset verbose = true showtiming = true ":min_tf :mayer" begin
         sol = solve(ocp, print_level = 0, tol = 1e-12)
-        @test sol.objective ≈ 2.0 rtol = 1e-2
+        @test objective(sol) ≈ 2.0 rtol = 1e-2
     end
 
     # min tf (lagrange)
@@ -34,7 +34,7 @@ function test_objective()
 
     @testset verbose = true showtiming = true ":min_tf :lagrange" begin
         sol = solve(ocp, print_level = 0, tol = 1e-12)
-        @test sol.objective ≈ 2.0 rtol = 1e-2
+        @test objective(sol) ≈ 2.0 rtol = 1e-2
     end
 
     # max t0 (free t0 and tf)
@@ -53,16 +53,16 @@ function test_objective()
 
     @testset verbose = true showtiming = true ":max_t0" begin
         sol = solve(ocp, print_level = 0, tol = 1e-12)
-        @test sol.objective ≈ 8.0 rtol = 1e-2
+        @test objective(sol) ≈ 8.0 rtol = 1e-2
     end
 
     @testset verbose = true showtiming = true ":max_t0 :explicit_grid" begin
         sol = solve(ocp, time_grid = LinRange(0, 1, 101), print_level = 0, tol = 1e-12)
-        @test sol.objective ≈ 8.0 rtol = 1e-2
+        @test objective(sol) ≈ 8.0 rtol = 1e-2
     end
 
     @testset verbose = true showtiming = true ":max_t0 :non_uniform_grid" begin
         sol = solve(ocp, time_grid = [0, 0.1, 0.6, 0.95, 1], print_level = 0, tol = 1e-12)
-        @test sol.objective ≈ 7.48 rtol = 1e-2
+        @test objective(sol) ≈ 7.48 rtol = 1e-2
     end
 end