You can use discretization of the problem into $N$ points, such that you only have to determine a finite number of parameters (assuming $f$ and $g$ are somewhat continuous functions). For the derivative and integration you can use Euler method, higher order methods can be used, but make the problem harder to solve.
The reformulation gives:
$$
h = \frac{t_1}{N-1}, \quad \vec{x} = [x_1, x_2, \dots, x_N], \quad \vec{y} = [y_1, y_2, \dots, y_N],
$$
$$
\begin{align}
\max_{\vec{x}, \vec{y}} & & \sum_{n=1}^{N-1} f(h (n-1), x_n, y_n) h\\
s.t. & & x_{n+1} = x_n + g(h (n-1), x_n, y_n) h, & & n = 1, 2, \dots, N-1
\end{align}
$$
You also have to add the boundary constraints to the equality constraints of the optimization problem. You can use multiple different methods to solve this problem, for example if you have access to Matlab you could use fmincon, which minimizes the cost function which can be fixed by adding an minus sign in front of the sum. Often you also have to supply an initial guess, which might also affect the solution, since different guesses might converge to different local maxima. By increasing $N$ you should get a more and more accurate solution, but it will probably take longer to solve. It might converge faster if you use the solution of a problem with less points and interpolate them and then use that as an initial guess for the problem of the larger number of points.