How can I have MATLAB solve a boundary value problem with free endpoints?

Suppose, I have the following boundary value problem:

f(x, xp, x2p, ..., xNp) = 0, subject to

g1(x(0), xp(0), ... x(N-1)p(0)) = 0. g2(x(tf), xp(tf), ... x(N-1)p(tf)) = 0.

where x?p denotes the ?th time derivative of a curve x: I --> M, where, M is a Riemannian manifold (You can assume M = R^n, for some natural number n for the purposes of this question.) The boundary conditions are given at the initial time by the function g1: I --> R^k, and the boundary conditions at the final time by the function g2: I --> R^m.

Now this problem can be solved by using the MATLAB function "bvp4c" as long as the final time tf is specified. The problem I would like to solve is when tf is free to vary. I will give an example of such a problem. Suppose we would like to intercept a curve with a particle in R^2. We still have a given differential equation f(x1,x2,x1p,x2p,...) = 0 to be satisfied, but we do not care where we intercept the path exactly, as long as it minimizes some functional. In this case, one of the coordinate functions of g2 would read something like

x(tf) = zeta(tf)

where zeta denotes the path we would like to intercept.

I believe I can solve this problem by an iterative application of "bvp4c", where I would guess a final time and solve the question with fixed endpoints. Afterwards, depending on the gradient of the unsatisfied boundary conditions, I can do a gradient based search to step the final time tf along the direction that reduces the norm of the unsatisfied boundary condition.

Is there an easier way of doing this? Any MATLAB functions anyone can suggest?

Best regards, Aykut