You can evaluate this using the Björck-Pereyra algorithm for solving Vandermonde systems, because you are evaluating $b^\top V^{-1}$ with $b=(2,0,\frac23,0,\frac25,0,\ldots)$, and the algorithm is known to be forward-stable (see Error analysis of the Björck-Pereyra algorithms for solving Vandermonde systems by Nick Higham, http://www.maths.manchester.ac.uk/~higham/narep/narep108.pdf).
Note: it seems that this analysis relies on the property that $0\leq x_1<x_2<\cdots <x_n$ (which is equivalent to $V$ being totally positive), as well as the elements of $b$ having alternating signs (which ensures there's no catastrophic cancellation in the subtractions below), in which case the errors are independent of the condition number, and it will still work in the more general case, without $0\leq x_1$, but the error bounds will be different. In any case, it takes $O(n^2)$ time, and avoids the problem with evaluating/integrating $L_i$'s, so it might be worth it even then, but I hadn't realized this point when I started writing this answer. You might be able to just map $x\mapsto \frac12(x+1)$ if that works for your problem.
I wrote a small Julia program to check that this actually works, and it gives $O(\epsilon_{\mathrm{mach}})$ relative errors.
module VandermondeInverse
using SpecialMatrices
function main(n=8)
X = Rational{BigInt}[k//(n-1) for k=0:n-1]
# X = convert(Vector{Rational{BigInt}}, linspace(-1, 1, n))
x = convert(Vector{Float64}, X)
A = convert(Matrix{Rational{BigInt}}, Vandermonde(X))
b = [i%2==0 ? 2//(i+1) : 0 for i=0:n-1]
println("Norm: ", norm(A, Inf))
println("Norm of inverse: ", norm(inv(A), Inf))
println("Condition number: ", cond(convert(Matrix{Float64}, A)))
ans = A'\b
println("True answer: ", ans)
B = convert(Matrix{Float64}, A)
c = convert(Vector{Float64}, b)
println("Linear solve: ", norm((B'\c - ans)./ans, Inf))
d = vec(c')
for k=1:n, l=n:-1:k+1
d[l] -= x[k]*d[l-1]
end
for k=n-1:-1:1, l=k:n
if l > k
d[l] /= x[l]-x[l-k]
end
if l < n
d[l] -= d[l+1]/(x[l+1] - x[l-k+1])
end
end
println("Neville elimination: ", norm((d-ans)./ans, Inf))
nothing
end
end
V = VandermondeInverse
Output:
julia> V.main(14)
Norm: 14.0
Norm of inverse: 1.4285962612120493e10
Condition number: 5.2214922998851654e10
True answer: Rational{Int64}[3202439130233//2916000,-688553801328731//52390800,19139253128382829//261954000,-196146528919726853//785862000,6800579086408939//11642400,-43149880138884259//43659000,32567483200938127//26195400,-7339312362348889//6237000,48767438804485271//58212000,-69618881108680969//157172400,44275410625421677//261954000,-2308743351566483//52390800,11057243346333379//1571724000,-209920276397//404250]
Linear solve: 1.5714609387747318e-8
Neville elimination: 1.3238218572356314e-15
If X isn't positive like in this test, then it seems the relative errors are of the same order as with a regular linear solve.
Why $b^\top V^{-1}$? It's actually a very useful common trick for working with polynomials of all types, but especially the Lagrange interpolating polynomials, converting the problem to a matrix form. The condition that defines $L_i$ is that $L_i(x_j)=\delta_{ij}$. Let $\alpha_{jk}$ be the coefficients of $L_k$, i.e.,
$$ L_k(x) = \sum_{j,k} \alpha_{j,k}x^j = (1,x,x^2,\ldots,x^n)^\top (\alpha_{0k},\ldots,\alpha_{nk}),
$$
and $L$ be the whole matrix of coefficients, arranged by columns:
$$ L = \begin{pmatrix} \alpha_{00}& \cdots & \alpha_{0n}\\\vdots &&\vdots\\ \alpha_{n0}&\cdots&\alpha_{nn} \end{pmatrix}. $$
Because of the definition of $L_k$ above as a vector product, multiplying $L$ on the left by $(1,x,\ldots,x^n)$ yields
$$ (1,x,x^2,\ldots,x^n)L = (L_0(x),L_1(x),\ldots,L_n(x)). $$
Using the condition $L_k(x_j)=\delta_{jk}$, this means that
$$ \begin{pmatrix} 1&x_0&x_0^2&\cdots&x_0^n\\ \vdots\\ 1&x_n&x_n^2&\cdots&x_n^n \end{pmatrix} L = I, $$
so $L = V^{-1}$, where $V$ is the Vandermonde matrix of $x_j$'s.
Finally, since $\int_{-1}^{1}x^k\,\mathrm{d}x = \frac{1+(-1)^k}{k+1}$, we have
$$ \int_{-1}^{1}L_k(x)\,\mathrm{d}x = \sum_j \alpha_{jk}\frac{1+(-1)^k}{k+1} = (2,0,\tfrac23,0,\tfrac25,0,\ldots)^\top (\alpha_{0k},\ldots,\alpha_{nk}).$$
So the $n+1$ numbers you are looking for, for $k=0\ldots n$, are given by $(2,0,\tfrac23,0,\ldots)^\top L$, where $L=V^{-1}$ is the inverse of the Vandermonde matrix.
