# Differential equation: mixed boundary condition, how to solve numerically?

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.

Instead of using optimisation, you can just solve a nonlinear equation. For each $$\alpha$$, evolve $$x_{\alpha},y_{\alpha}$$ by solving the corresponding differential equation (with initial condition $$y(0)=\alpha$$), and define $$F(\alpha) = y_{\alpha}(1) - g(x_{\alpha}(1))$$. Now you "just" need to solve $$F(\alpha)=0$$. I believe you can find inspiration in the "shooting method".
I believe that you have no much better choice than trial & error, but this can be recast as a standard equation resolution problem: the unknown is $$y(0)$$ and the target function is $$y(1)-g(x(1))$$, where $$x(1)$$ and $$y(1)$$ are obtained by numerical integration, starting from $$x(0),y(0)$$.