I will make the following assumptions

- you can evaluate your function wherever your like.
- you are unconcerned with machine precision (this should not be an issue for the number of coefficients you are after).

Then the solution is to use an n-point Gauss-Legendre quadrature in order to optimally extract n Legendre coefficients.

To obtain this result, the key point is that the function space spanned by the first $n$ Legendre polynomials is simply the space of the zero-th to n-th polynomial moments. More specifically:

$$ \mathbb{P}_n\equiv\mathrm{span}(P_0(x),P_1(x),P_2(x),\ldots,P_n(x)) \equiv \mathrm{span}(1,x,x^2,\ldots,x^n).$$

Hence to perfectly integrate the $\{0,\ldots,n-1\}$ Legendre moments, we are really looking to integrate an order-$2n-2$ polynomial exactly. We are looking for an answer to the following.

**Question.** Given that $f$ is a polynomial of at most order $n$, is there a numerical integration scheme that will integrate the product $f$ exactly?

In this form, the question has a known, provable solution.

**Theorem 1 (Gauss-Legendre Quadrature).** The $n$-node Gaussian quadrature scheme, whose nodes are defined at the roots of the order $n$ Legendre polynomial, and whose weights are defined via Lagrange interpolants, will peform the integration exactly$$\int_0^1f(x)\,dx=\sum_{i=1}^m w_i f(x_i)$$for all $f\in\mathbb{P}_{2n-1}$, where $\mathbb{P}_{2n-1}$ denotes the space of order $2n-1$ polynomials.

In terms of where to find code for Gauss-Legendre Quadrature, I would look for it on MATLAB central file exchange.