In the specific problem you ask about (unlike Richard's more general answer), it turns out that you can relax the $=$ constraint into a $\leq$ without changing the optimal value, and that the resulting convex problem can be solved with the CVX software Richard mentioned.
Details: Intuitively the relaxation is possible since if the function had arc length strictly less than $L$, you could scale it up to have arc length $L$ exactly, but scaling the function up will necessarily increase the integral. As a result, if we solve the optimization problem with a $\leq$ replacing the $=$ in the constraint, the optimal value and function $y(x)$ is the same.
It turns out the resulting problem is convex in the parameter $y(x)$, albeit infinite dimensional. The objective is linear. Moreover, the function $\sqrt{1 + y'(z)^2} = \|(1, y'(z))\|_2$ is a convex function of $y(x)$ as a norm of a linear function of $y(x)$ (the derivative is linear). Thus, since the integral of a family of convex functionals is convex, the mapping $y(x) \mapsto \int_a^b \sqrt{1 + y'(x)^2}\,dx$ is convex. The constraint is then convex, as it asks for $y(x)$ to lie in a sub-level set of a convex function of $y(x)$.
The following Python code gives an example of solving this problem with the CVXPY software Richard mentioned:
import numpy as np
import cvxpy as cp
# Problem Data
n = 101 # points in interval [0,1] to discretize y at
L = 1.5
y = cp.Variable(n) # y[i] = y(i/n) for i=0,1,...,n
x = np.arange(n) / n
dy = y[1:] - y[:-1]
dx = x[1:] - x[:-1]
# Construct Optimization Problem
objective = (cp.sum(y[:-1]) + cp.sum(y[1:])) / 2 / n
constraints = [
y[0] == 0, y[-1] == 0,
cp.sum(cp.norm2(cp.vstack([dy, dx]), axis=0)) <= L,
]
prob = cp.Problem(cp.Maximize(objective), constraints)
# Solve Problem
prob.solve(verbose=True)
The solution (found at y.value) is given as follows:
As Richard points out, you can't hope to be able to add generic constraints to this problem and continue to have access to a simple and/or efficient solver. Nevertheless, you can add convex constraints like this arc-length one in a pretty straightforward way without too much trouble. For example, the following gives the output of the same optimization problem if I now add the constraint $y(1/2) \geq 0.55$:
