Symbolic solution for the energy of potential flow

In a numerical model, energy is conserved with some accuracy; in this example, the deviation from the initial value is about 1.5*10^-12

m = 1; omega0 = 1; eps = 1/100; v0 = 1;
eq = m*x''[t] == -m*(omega0^2*x[t] + eps*x[t]^3);
ic = {x[0] == 0, x'[0] == v0};

X = NDSolveValue[{eq, ic}, x, {t, 0, 10}, WorkingPrecision -> 30];

Plot[m/2*X'[t]^2 + m/2*omega0^2*X[t]^2 + m/4*eps*X[t]^4 - 
  m/2*v0^2, {t, 0, 10},AxesLabel -> {"t", "E-E0"}]

Using the law of conservation of energy, we express $x'(t) $ and then time as a function of $x$

t=Integrate[1/Sqrt[v0^2 - omega0^2*x^2 - eps/2*x^4], x]
(*-((I Sqrt[2 + (2 eps x^2)/(omega0^2 - Sqrt[omega0^4 + 2 eps v0^2])]
   Sqrt[1 + (eps x^2)/(omega0^2 + Sqrt[omega0^4 + 2 eps v0^2])]
   EllipticF[
   I ArcSinh[Sqrt[eps/(omega0^2 - Sqrt[omega0^4 + 2 eps v0^2])] x], (
   omega0^2 - Sqrt[omega0^4 + 2 eps v0^2])/(
   omega0^2 + Sqrt[omega0^4 + 2 eps v0^2])])/(
 Sqrt[eps/(omega0^2 - Sqrt[omega0^4 + 2 eps v0^2])] Sqrt[
  2 v0^2 - 2 omega0^2 x^2 - eps x^4]))*)

This problem can be solved symbolically as follows. Multiply the expression (m (omega0^2 x[t] + eps x[t]^3) + m x''[t]) by x'[t] and integrate to obtain an expression for the energy of this nonlinear oscillator.

eq = Integrate[(m (omega0^2 x[t] + eps x[t]^3) + m x''[t]) x'[t], t]
(* 1/2 m omega0^2 x[t]^2 + 1/4 eps m x[t]^4 + 1/2 m x'[t]^2 *)

with constant of integration v0, the conserved energy. Then, apply DSolve.

s = DSolve[eq == v0, x[t], t] // Last
(* {x[t] -> InverseFunction[-((I EllipticF[I ArcSinh[Sqrt[(eps m)/(m omega0^2 -
    Sqrt[m (m omega0^4 + 4 eps v0)])] #1], (m omega0^2 - Sqrt[m (m omega0^4 + 
    4 eps v0)])/(m omega0^2 + Sqrt[m (m omega0^4 + 4 eps v0)])] Sqrt[1 + (eps m #1^2)
    /(m omega0^2 - Sqrt[m (m omega0^4 + 4 eps v0)])] Sqrt[1 + (eps m #1^2)
    /(m omega0^2 + Sqrt[m (m omega0^4 + 4 eps v0)])])/(Sqrt[(eps m)/(m omega0^2 - 
    Sqrt[m (m omega0^4 + 4 eps v0)])] Sqrt[4 v0 - m #1^2 (2 omega0^2 + eps #1^2)])) &]
    [t/(Sqrt[2] Sqrt[m]) + C[1]]} *)

(The other solution is the negative of the first.) Since the question requests t as a function of x, s must be inverted. In the absence of a Mathematica command to accomplish this, we use the following ungainly expression.

st = Rule[(s[[1, 2, 1]] /. C[1] -> 0) Sqrt[2] Sqrt[m], 
    Head[s[[1, 2]]][[1]][x[t]] Sqrt[2] Sqrt[m]]
(* t -> -((I Sqrt[2] Sqrt[m] EllipticF[I ArcSinh[
       Sqrt[(eps m)/(m omega0^2 - Sqrt[m (m omega0^4 + 4 eps v0)])] x[t]], 
       (m omega0^2 - Sqrt[m (m omega0^4 + 4 eps v0)])/
       (m omega0^2 + Sqrt[m (m omega0^4 + 4 eps v0)])] 
       Sqrt[1 + (eps m x[t]^2)/(m omega0^2 - Sqrt[m (m omega0^4 + 4 eps v0)])] 
       Sqrt[1 + (eps m x[t]^2)/(m omega0^2 + Sqrt[m (m omega0^4 + 4 eps v0)])])/(
       Sqrt[(eps m)/(m omega0^2 - Sqrt[m (m omega0^4 + 4 eps v0)])] 
       Sqrt[4 v0 - m x[t]^2 (2 omega0^2 + eps x[t]^2)])) *)  

This result for various values of eps can be plotted as

st /. {m -> 1, omega0 -> 1, v0 -> 1};
Plot[Evaluate@Table[Last[%], {eps, {10^1, 1, 10^-1, 10^-10}}], {x[t], -2, 2}, 
    AxesLabel -> {x, t}, AspectRatio -> 1, ImageSize -> Large, 
    LabelStyle -> {Bold, Black, 15}]

Decreasing eps corresponds to increasing values of x and t at the turning points.