diff --git a/.readthedocs.yml b/.readthedocs.yml
index f9514ee9..d514e60c 100644
--- a/.readthedocs.yml
+++ b/.readthedocs.yml
@@ -1,8 +1,8 @@
 version: 2
 build:
-  os: "ubuntu-22.04"
+  os: "ubuntu-24.04"
   tools:
-    python: "mambaforge-22.9"
+    python: "miniconda-latest"
 sphinx:
   configuration: docs/conf.py
 formats: all
diff --git a/README.rst b/README.rst
index 54ad8c1e..235221c2 100644
--- a/README.rst
+++ b/README.rst
@@ -79,6 +79,7 @@ The required dependencies are as follows:
 - numpy >= 1.19.0
 - python 3.9-3.13
 - scipy >= 1.5.0
+- setuptools
 - sympy >= 1.6.0
 
 The optional dependencies are as follows:
@@ -136,7 +137,7 @@ set the ``LD_LIBRARY_PATH`` so that you can link to Ipopt when installing
 
 Once Ipopt is installed and accessible, install conda then create an environment::
 
-   $ conda create -n opty-custom -c conda-forge cython numpy pip scipy sympy
+   $ conda create -n opty-custom -c conda-forge cython numpy pip scipy setuptools sympy
    $ source activate opty-custom
    (opty-custom)$ pip install cyipopt  # this will compile cyipopt against the available ipopt
    (opty-custom)$ pip install opty
@@ -145,7 +146,7 @@ If you want to develop opty, create a conda environment with all of the
 dependencies installed::
 
    $ conda config --add channels conda-forge
-   $ conda create -n opty-dev python sympy numpy scipy cython ipopt cyipopt matplotlib pytables pydy pandas pytest sphinx sphinx-gallery numpydoc
+   $ conda create -n opty-dev python sympy numpy scipy cython ipopt cyipopt matplotlib pytables pydy pandas pytest setuptools sphinx sphinx-gallery numpydoc
    $ source activate opty-dev
 
 Next download the opty source files and install with::
diff --git a/opty-dev-env.yml b/opty-dev-env.yml
index e2e073a0..e995ba18 100644
--- a/opty-dev-env.yml
+++ b/opty-dev-env.yml
@@ -13,6 +13,7 @@ dependencies:
   - pytest-cov
   - python
   - scipy >=1.5.0
+  - setuptools
   - sphinx
   - sphinx-gallery
   - sympy >=1.13  # version required for the docs
diff --git a/opty/direct_collocation.py b/opty/direct_collocation.py
index c598695e..47059a0f 100644
--- a/opty/direct_collocation.py
+++ b/opty/direct_collocation.py
@@ -128,6 +128,8 @@ class Problem(cyipopt.Problem):
                       eomn2, ... eomnN,
                       c1, ..., co]
 
+    The attributes may be accessed as follows: ``Problem_instance.collocator.name_of_attribute``
+
     """
 
     INF = 10e19
@@ -340,7 +342,7 @@ def jacobian(self, free):
 
         Returns
         =======
-        jac_vals : ndarray, shape((2*n + q + r + s)*(n*(N - 1) + o), )
+        jac_vals : ndarray, shape((2*n + q + r + s)*(n*(N - 1)) + o, )
             Non-zero Jacobian values in triplet format.
 
         """
@@ -597,6 +599,44 @@ def plot_objective_value(self):
 
         return ax
 
+    def parse_free(self, free):
+        """Parses the free parameters vector and returns it's components.
+
+        Parameters
+        ==========
+        free : ndarray, shape(n*N + q*N + r + s)
+            The free parameters of the system.
+
+        Returns
+        =======
+        states : ndarray, shape(n, N)
+            The array of n states through N time steps.
+        specified_values : ndarray, shape(q, N) or shape(N,), or None
+            The array of q specified inputs through N time steps.
+        constant_values : ndarray, shape(r,)
+            The array of r constants.
+        time_interval : float
+            The time between collocation nodes. Only returned if
+            ``variable_duration`` is ``True``.
+
+        Notes
+        =====
+
+        - N : number of collocation nodes
+        - n : number of unknown state trajectories
+        - q : number of unknown input trajectories
+        - r : number of unknown parameters
+        - s : number of unknown time intervals (s=1 if ``variable duration`` is
+          ``True`` else s=0)
+
+        """
+
+        n = self.collocator.num_states
+        N = self.collocator.num_collocation_nodes
+        q = self.collocator.num_unknown_input_trajectories
+        variable_duration = self.collocator._variable_duration
+
+        return parse_free(free, n, q, N, variable_duration)
 
 class ConstraintCollocator(object):
     """This class is responsible for generating the constraint function and the
@@ -604,6 +644,81 @@ class ConstraintCollocator(object):
     for a non-linear programming problem where the essential constraints are
     defined from the equations of motion of the system.
 
+    Attributes
+    ==========
+    current_discrete_state_symbols : n-tuple
+        The symbols for the current discrete states.
+    current_discrete_specified_symbols : q-tuple
+        The symbols for the current discrete specified inputs.
+    discrete_eom : sympy.Matrix, shape(n, 1)
+        Discretized equations of motion. Depending on the integration method
+        used.
+    eom: sympy.Matrix, shape(n, 1)
+        The equations of motion used.
+    input_trajectories : tuple
+        known_input_trajectories + unknown_input_trajectories.
+    instance_constraints : o-tuple
+        The instance constraints used in the optimization.
+    integration_method : str
+        The integration method used.
+    known_parameters : tuple
+        The symbols of the known parameters in the problem.
+    known_parameter_map : dict
+        A mapping of known parameters to their values.
+    known_trajectory_map : dict
+        A mapping of known trajectories to their values.
+    known_trajectory_symbols : (m-q)-tuple
+        The known trajectory symbols.
+    next_discrete_specified_symbols : q-tuple
+        The symbols for the next discrete specified inputs.
+    next_discrete_state_symbols : n-tuple
+        The symbols for the next discrete states.
+    node_time_interval : float or sympy.Symbol
+        The time interval between the collocation nodes. float if the interval
+        is fixed, ``sympy.Symbol`` if the interval is variable.
+    num_collocation_nodes : int
+        Number of times spaced evenly between the initial and final time of
+        the optimization = N.
+    num_constraints : int
+        The number of constraints = (num_collection_nodes-1)*num_states +
+        len(instance_constraints).
+    num_free : int
+        Number of variables to be optimized = n*N + q*N + r + s.
+    num_input_trajectories : int
+        The number of input trajectories = len(input_trajectories).
+    num_instance_constraints : int
+        The number of instance constraints = len(instance_constraints).
+    num_known_trajectories : int
+        The number of known trajectories = len(known_trajectory_symbols).
+    num_parameters : int
+        The number of parameters = len(parameters).
+    num_states : int
+        The number of states = len(state_symbols) = n.
+    num_unknown_input_trajectories : int
+        The number of unknown input trajectories =
+        len(unknown_input_trajectories).
+    num_unknown_parameters : int
+        The number of unknown parameters = r.
+    parameters : tuple
+        known_parameters + unknown_parameters.
+    parallel : bool
+        Whether to use parallel processing or not.
+    previous_discrete_state_symbols : n-tuple
+        The symbols for the previous discrete states.
+    show_compile_output : bool
+        Whether to show the compile output or not.
+    state_derivative_symbols : n-tuple
+        symbols for the time derivatives of the states.
+    time_symbol : sympy.Symbol
+        The symbol used to represent time, usually `t`.
+    tmp_dir
+        The temporary directory used to store files generated.
+    unknown_input_trajectories : q-tuple
+        The unknown input trajectories symbols.
+    unknown_parameters : r-tuple
+        The unknown parameters in the problem, in the sequence in which they
+        appear in the solution of the optimization.
+
     Notes
     =====
 
@@ -618,8 +733,12 @@ class ConstraintCollocator(object):
     - nN + qN + r + s : number of free variables
     - n(N - 1) + o : number of constraints
 
-    """
+    Some of the attributes are explained in more detail under Parameters below.
+
+    It is best to treat ``ConstraintCollocator`` as immutable, changing
+    attributes after initialization will inevitably fail.
 
+    """
     def __init__(self, equations_of_motion, state_symbols,
                  num_collocation_nodes, node_time_interval,
                  known_parameter_map={}, known_trajectory_map={},
diff --git a/opty/tests/test_direct_collocation.py b/opty/tests/test_direct_collocation.py
index adc2f3f8..ee598876 100644
--- a/opty/tests/test_direct_collocation.py
+++ b/opty/tests/test_direct_collocation.py
@@ -10,7 +10,7 @@
 from pytest import raises
 
 from ..direct_collocation import Problem, ConstraintCollocator
-from ..utils import create_objective_function, sort_sympy
+from ..utils import create_objective_function, sort_sympy, parse_free
 
 
 def test_pendulum():
@@ -1566,3 +1566,124 @@ def test_for_algebraic_eoms():
         )
 
     assert excinfo.type is ValueError
+
+
+def test_prob_parse_free():
+    """
+    Test for parse_free method of Problem class.
+    ===========================================
+
+    This test whether the parse_free method of the Problem class works as
+    the parse_free in utils.
+
+    **States**
+
+    - :math:`x_1, x_2, ux_1, ux_2` : state variables
+
+    **Control**
+
+    - :math:`u_1, u_2` : control variable
+
+    """
+
+    t = mech.dynamicsymbols._t
+
+    x1, x2, ux1, ux2 = mech.dynamicsymbols('x1 x2 ux1 ux2')
+    u1, u2 = mech.dynamicsymbols('u1 u2')
+    h = sym.symbols('h')
+    a, b = sym.symbols('a b')
+
+    # equations of motion.
+    # (No meaning, just for testing)
+    eom = sym.Matrix([
+            -x1.diff(t) + ux1,
+            -x2.diff(t) + ux2,
+            -ux1.diff(t) + a*u1,
+            -ux2.diff(t) + b*u2,
+    ])
+
+    # Set up and Solve the Optimization Problem
+    num_nodes = 11
+    t0, tf = 0.0, 0.9
+    state_symbols = (x1, x2, ux1, ux2)
+    control_symbols = (u1, u2)
+
+    interval_value = (tf - t0)/(num_nodes - 1)
+    times = np.linspace(t0, tf, num_nodes)
+
+    # Specify the symbolic instance constraints, as per the example.
+    instance_constraints = (
+        x1.func(t0) - 1.0,
+        x2.func(t0) + 1.0,
+    )
+
+    # Specify the objective function and form the gradient.
+
+    def obj(free):
+        return sum([free[i]**2 for i in range(2*num_nodes)])
+
+    def obj_grad(free):
+        grad = np.zeros_like(free)
+        grad[:2*num_nodes] = 2*free[:2*num_nodes]
+        return grad
+
+    # Create the optimization problem and set any options.
+    prob = Problem(
+            obj,
+            obj_grad,
+            eom,
+            state_symbols,
+            num_nodes,
+            interval_value,
+            instance_constraints=instance_constraints,
+)
+
+    # Give some estimates for the trajectory.
+    initial_guess = np.random.rand(prob.num_free)
+    initial_guess1 = initial_guess
+
+    # check whether same results.
+    statesu, controlsu, constantsu = parse_free(initial_guess1,
+            len(state_symbols), len(control_symbols), num_nodes)
+
+    states, controls, constants = prob.parse_free(initial_guess)
+    np.testing.assert_allclose(states, statesu)
+    np.testing.assert_allclose(controls, controlsu)
+    np.testing.assert_allclose(constants, constantsu)
+
+    # test with variable interval_value
+    interval_value = h
+    t0, tf = 0.0, (num_nodes - 1)*interval_value
+    def obj(free):
+        return sum([free[i]**2 for i in range(2*num_nodes)])
+
+    def obj_grad(free):
+        grad = np.zeros_like(free)
+        grad[:2*num_nodes] = 2*free[:2*num_nodes]
+        return grad
+
+    # Create the optimization problem and set any options.
+    prob = Problem(
+            obj,
+            obj_grad,
+            eom,
+            state_symbols,
+            num_nodes,
+            interval_value,
+            instance_constraints=instance_constraints,
+)
+
+    # Give some estimates for the trajectory.
+    initial_guess = np.random.rand(prob.num_free)
+    initial_guess1 = initial_guess
+
+    # check whether same results.
+    statesu, controlsu, constantsu, timeu = parse_free(initial_guess1,
+        len(state_symbols), len(control_symbols),
+        num_nodes, variable_duration=True)
+
+    states, controls, constants, times = prob.parse_free(initial_guess)
+    np.testing.assert_allclose(states, statesu)
+    np.testing.assert_allclose(controls, controlsu)
+    np.testing.assert_allclose(constants, constantsu)
+    np.testing.assert_allclose(timeu, times)
diff --git a/setup.py b/setup.py
index 34f7b404..33d3eefc 100644
--- a/setup.py
+++ b/setup.py
@@ -20,6 +20,7 @@
         'cython>=0.29.19',
         'numpy>=1.19.0',
         'scipy>=1.5.0',
+        'setuptools',  # provides distutils for Python >=3.13
         'sympy>=1.6.0',
     ],
     extras_require={