I have a differential equation (dot means derivative w.r.t time) $$(\dot x, \dot y) = f(x,y)$$ and I am given the initial condition for $x$, but a final condition for $y$: $$x(0),\qquad y(1)=g(x(1))$$ The second simply means that the final value of $y$ depends on the final value of $x$ with no known initial value.
Since I do not have an initial condition for $y$, my very crude attempt to solve it numerically is to pick a random $y(0)$, evolve the system in time, and compare the final value of $y$, $y(1)$, with the value it should have:
$$ I(y(0)) = \vert y(1)-g(x(1)) \vert ^2 $$
So when $y(1)=g(x(1))$ I will fulfill my condition. $I$ would be a non-negative "cost" associated with each initial condition. I applied a gradient descent numerically but this is quite unstable, meaning it doesn't converge. I also applied black-box optimizers based on Newton's method (the Broyden–Fletcher–Goldfarb–Shanno algorithm) but it very rarely converges.
Are there any other algorithms for solving mixed boundary conditions of this type?
I'm also open to "rules of thumb" which, while not justified mathematically, can make solving optimization problems easier.
Context: I am trying to solve an optimal control problem. I apply Pontryagin's theorem and I get the mixed boundary conditions. $x$ and $y$ are actually vectors of dimension 4, but the above example captures the difficulty well enough in 2D.