I define this Lagrangian describing a physical system under investigation (two pendula connected by non-linear springs):

    Lag = 1/2 m (x1'[t]^2 + x2'[t]^2) - (m g)/(
    2 l) (x1[t]^2 + x2[t]^2) - (k/
      2 (x1[t]^2 + (x2[t] - x1[t])^2 + (x2[t] - x1[t])^2 + (L - 
      x2[t])^2) + \[Alpha] (x1[
     t]^3 + (x2[t] - x1[t])^3 + (x2[t] - x1[t])^3 + (L - 
      x2[t])^3) + \[Beta] (x1[
     t]^4 + (x2[t] - x1[t])^4 + (x2[t] - x1[t])^4 + (L - 
      x2[t])^4));

where all coefficients are constant. I ask Mathematica to calculate the Euler equations as follows

    ee = EulerEquations[Lag, {x1[t], x2[t]}, t]

Hence I ask M to solve the corresponding non linear system of two second order ODEs, using NDSolve

    sol = NDSolve[
    Join[ee, {x1[0] == 0, x2[0] == 0, x1'[0] == 0.1, 
    x2'[0] == 0}], {x1[t], x2[t]}, {t, 0, 10000}][[1, All, 2]]

NDSolve works pretty fast. But after that I need to calculate the average energy of each body as a function of time and I use those two list of lists:

    u1=Table[NIntegrate[(m g)/(2 l) sol[[1]]^2 + 
    k/2 (sol[[1]]^2 + (sol[[2]] - sol[[1]])^2) + \[Alpha] (sol[[
    1]]^3 + (sol[[2]] - sol[[1]])^3) + \[Beta] (sol[[
    1]]^4 + (sol[[2]] - sol[[1]])^4), {t, 0, i}], {i, n}]];
    uMedia1=Accumulate[u1];

It takes too much to plot sol[[1]] (or 2) and even longer to NIntegrate over each time interval (by the way Accumulate is very fast). As a consequence I guess I'm doing something wrong with that.

Is this the correct procedure to solve this problem? I appreciate your help. Marta

k = 1; l = 1; L = 2; m = 1; g = 9.8; [Alpha] = 0.5; [Beta] = 0;
[Lambda] = 0.1 ;