Thanks to Dimitry, I have edited the question to correct some errors.
I have been working on the optimal control of an underwater vehicle attempting to dock onto a charging station. Limiting the problem to three degrees of freedom, the nonlinear dynamic equations can be expressed as
$\mathbf{\dot{x}}(t) = \mathbf{g}(\mathbf{x}(t),\mathbf{u}(t)),$
where $\mathbf{u} = \begin{bmatrix} Q & \delta_\mathrm{s} \end{bmatrix}^T$ and $\mathbf{x} = \begin{bmatrix} x & z & \theta & u & w & q & n & u_\mathrm{p} \end{bmatrix}^T$.
In order to dock, the state vector should have the following initial and end conditions
$\mathbf{x}_0=\begin{bmatrix} x_0 & z_0 & \theta_0 & u_0 & w_0 & q_0 & n_0 & u_{\mathrm{p},0} \end{bmatrix}^T,$
$\mathbf{x}_\mathrm{f}=\begin{bmatrix} x_\mathrm{f} & z_\mathrm{f} & \theta_\mathrm{f} & u_\mathrm{f} & w_\mathrm{f} & q_\mathrm{f} & \mathrm{free} & \mathrm{free} \end{bmatrix}^T.$
Note that the last two states are free at the end time $t_\mathrm{f}$ and I am not sure if this will cause issues.
Taking into account that high values of the control input are undesirable, I have come up with the following cost-function for the constrained optimisation problem (due to physical limits):
$J = 0.5 \int_0^{t_\mathrm{f}} \sum_{j=1}^2 w_{u,j} \left[u_j(t) \right]^2 \mathrm{d}t $ subject to
$\mathbf{\dot{x}}(t) = \mathbf{g}(\mathbf{x}(t),\mathbf{u}(t)),$
$\mathbf{x}(0) = \mathbf{x}_0,$
$\mathbf{x}(t_\mathrm{f}) = \mathbf{x}_\mathrm{f},$
$\mathbf{u}_\mathrm{min} \leq \mathbf{u}(t) \leq \mathbf{u}_\mathrm{max}.$
Note that at the moment, I am considering a fixed end time.
Using an indirect method, the costate equations may be found as follows.
Based on my cost function, the Hamiltonian is
$H(t) = -0.5\left( w_{u,1}Q^2(t) + w_{u,2}\delta_\mathrm{s}^2(t) \right) + \mathbf{p}^T(t) \mathbf{g}(\mathbf{x}(t),\mathbf{u}(t)),$
so that I get the following equations for the state and costate
$\mathbf{\dot{x}}(t) = \frac{\partial H(\mathbf{x},\mathbf{u},\mathbf{p},t)}{\partial \mathbf{p}} ,$
$\mathbf{\dot{p}}(t) = -\frac{\partial H(\mathbf{x},\mathbf{u},\mathbf{p},t)}{\partial \mathbf{x}} ,$
with the stationary condition
$\frac{\partial H(\mathbf{x},\mathbf{u},\mathbf{p},t)}{\partial \mathbf{u}} = 0.$
Furthermore, the boundary and transversality conditions can be expressed as
$\mathbf{x}(0) = \mathbf{x}_0, $
$p_n (t_\mathrm{f}) = 0,$
$p_{u_\mathrm{p}} (t_\mathrm{f}) = 0.$
My first question is how I can convert this into a boundary value problem so that I can use a tool like bvp4c in MATLAB to solve it. In particular, I know that the state equations should be solved forward in time and the costate backwards in time, but am not sure how to do that with bvp4c. My other problem is that I feel that I am not using the values of $\mathbf{x}_\mathrm{f}$ with fixed values. As a result, the number of boundary conditions is less than the number of ordinary differential equations, which is a problem.
Then, I have problems converting the problem into a free end time problem. I have seen that the procedure for this involves rescaling the time to $\tau=\frac{t}{t_\mathrm{f}}$, with $\tau \in [0,1]$. As a result, the optimal control problem is modified to
$J = 0.5 \int_0^{t_\mathrm{f}} \left\lbrace \sum_{j=1}^2 w_{u,j} \left[u_j(t) \right]^2 + 2 \right\rbrace \mathrm{d}t $ subject to
$\mathbf{\dot{x}}_\mathrm{new}(t) = \begin{bmatrix} \mathbf{\dot{x}}(t) \\ \dot{t} \end{bmatrix}= \begin{bmatrix} t_\mathrm{f} \mathbf{g}(\mathbf{x}(\tau),\mathbf{u}(\tau)) \\ 0 \end{bmatrix},$
$\begin{bmatrix}\mathbf{x}(0) \\ t(0) \end{bmatrix}= \begin{bmatrix}\mathbf{x}_0 \\ t_\mathrm{f} \end{bmatrix},$
$\begin{bmatrix}\mathbf{x}(1) \\ t(1) \end{bmatrix}= \begin{bmatrix}\mathbf{x}_\mathrm{f} \\ t_\mathrm{f} \end{bmatrix},$
$\mathbf{u}_\mathrm{min} \leq \mathbf{u}(\tau) \leq \mathbf{u}_\mathrm{max},$
$0<\tau<1.$
The new Hamiltonian then becomes
$H(t) = - 1 - 0.5\left( w_{u,1}Q^2(\tau) + w_{u,2}\delta_\mathrm{s}^2(\tau) \right) + \mathbf{p}^T(\tau) t_\mathrm{f} \mathbf{g}(\mathbf{x}(\tau),\mathbf{u}(\tau))$ and the stationarity conditions
$\mathbf{\dot{p}}(\tau) = -t_\mathrm{f} \frac{\partial H(\mathbf{x},\mathbf{u},\mathbf{p},\tau)}{\partial \mathbf{x}},$
including the additional costate for the time.
My second question is whether this procedure for free end time is now expressed error-free. Additionally, as in the first question, how can I express this as a boundary value problem I can solved with bvp4c?
Many thanks in advance for the help!