@Al86 

Look at this ansatz_2_mmcdara.mw

In addition: here is a A_Toy_Problem.mw which explains how the perturbation method operates and how to perform it step-by-step in Maple.

Note: I used the package LargeExpressions in this filr, not because it is necessary (you coud use coeff/coeffs instead), but to help you understand how I proceeeded in file ansatz_2_mmcdara.mw.

It seems method=gear is the method that best conserves the total energy  (absolute error less than 10^-9 over the range t=0..100, so about 10³ times smaller than Rosenbrock)

Hamiltoninan_procedure_mmcdara.mw

About Gear's method:
In French Hairer   (Gear's method is named BDF in this paper)
In English Novik et al. (specificaly chapter V. PERFORMANCE OF THE ALGORITHMS)

More generaly:
If you are interested in energy conserving methods for hamitonian systems I advice you to give a look to this introductory Wikipedia  article.

At last: your ode contains the term k*diff(x(t), t) that you cancel out setting k=0.
I guess you are aware that the ode no longer defines a conservative system as soon as k <> 0...

I'm not sure the ansatz you use is is the one that should be used.

In the attached file, up to relation (ç), you will found how to get an expression of the form

sum(A[k]*epsilon^k, k=M..N)

where each A[k] represents the kth order differential equation.
With the ansatz you chosed M=-9 and N=11.

From relation (10) to relation (19) the solutions x[-1](t), x[0](t), x[1](t) are computed.
You will see they all have the sape formal expression p+q*exp(-t), and thus is  x(t) writes p+q*exp(-t) too.

It's hard to go further without having the initial conditions for x[-1](t), x[0](t), x[1](t), but  it seems the ansatz you chosed is not well suited for the original equation.
Indeed, assuming x(t) <> 0 for each t in the integration domain, the ode can be transform into what I named ode_2.
Collecting ode_2  wrt epsilon writes 

sum(B[k]*epsilon^k, k=0..4)

So why did you used an ansatz containing x[-1](t)/epsilon?

ansatz_1_mmcdara.mw

@C_R 

and thank you for the comment.
Doing this kind animation was quite funny and I did a few others in a row. The first one was quite time-consuming, but once you have understood how plottools (mainly) works it is not that difficult to do amazing things.

I completely agree with your last sentence, even though I'm no longer involved in using MapleSim.


A story:
For years, I've been trying to promote MapleSim in my company. Until it made a deal with a big actor (Dassult Systèmes) which provides turnkey solutions, or so they say. In particular  a product named Dymola (based on Modelica) that some colleagues wrongly considered as an alternative to MapleSim (they just saw the pictures and movies were quite close but forgot the CAS aspect behind).
Whatever, when it goes to big money MapleSoft doesn't hold a candle to aggressive distributors, specially when outsourcing becomes the rule.
In the end MapleSim was rejected and we are now dependent on Dassault's goodwill to develop, or not, what we really need. End of the story.

@JAMET 

here3.mw

@Christian Wolinski 

First reference
Gauth  ( same problem and same result of yours [which is obviously correct])

Second reference
Doubtnut which gives answer 16
But this is not the same problem and @yangtheary made a mistake while reformulating the problem: the term power set makes all the difference!
(note the typo {x|x| < 3 ...} instead of  {x : |x| < 3 ...} on the screenshot below [the explanations the previous link provides prove t is indeed a typo]  }
 

 
In case @yangtheary would know how Maple can provide a solution, here is a way
 

Download With_Maple.mw

@JAMET 

write point(F, [2.329411765, -2.567510609]) if you want to check that AreHarmonic(Q, F, P, M) returns true.
Neither AreHarmonic nor AreCollinear support floats.

Redefining point F by 

X := 198/85:              # X := convert(2.329411765, rational);
Y := -(126/85)*sqrt(3):   # see the attached file
point(F, [X, Y]):

remove all floats in points Q, F, P, M and make them collinear.

Nevertheless (still look to the attached file), (Q, F, P, M) are not harmonic:

AreHarmonic(Q, F, P, M); 
                             false

# CrossRatio(Q, F, P, M) should be equal to -1

CrossRatio(Q, F, P, M):
evalf(%);
                          1.999999999

So (Q, F, P, M) are definitely not harmonic.

But...
 (Q, P, F, M) are  harmonic:

AreHarmonic(Q, P, F, M);
                             true
evalf[20](CrossRatio(Q, P, F, M));
                     -1.0000000000000000001

here2.mw

Solution here2.txt  (open the link and, if you want, copy-paste its content within a void Maple worksheet)

restart:
kernelopts(version);
OneFrame := proc(k)
local Courbe, T, a, b, c, t, P, Q, NormM, F, Ell, sol, N1, N2, dr, tx;
local MySol:
uses geometry:
_EnvHorizontalName := x: _EnvVerticalName := y:
a := 11; b := 7; c := sqrt(a^2 - b^2); t := 1/3*Pi;
Ell := x^2/a^2 + y^2/b^2 = 1;
point(T, (a^2 - b^2)*cos(t)^3/a, -(a^2 - b^2)*sin(t)^3/b);
Courbe := plots:-implicitplot(Ell, x = -a - 10 .. a + 10, y = -b - 10 .. b + 10, scaling = constrained, color = blue, grid=[50, 50]);
NormM := plots:-implicitplot(y - b*sin(t) = a*sin(t)*(x - a*cos(t))/(b*cos(t)), x = -a - 5 .. a + 10, y = -b - 10 .. b + 10, color = orange);
line(Per, y - b*sin(t) = a*sin(t)*(x - a*cos(t))/(b*cos(t)), [x, y]);
point(P, [solve(subs(y = 0, lhs(Equation(Per)))), 0]);
point(Q, [0, solve(subs(x = 0, lhs(Equation(Per))))]);
point(M, [a*cos(t), b*sin(t)]);
point(N1, [a*cos(k), b*sin(k)]);
point(F, [2.329411765, -2.567510609]);
line(L, [N1, F]);
sol := solve({Equation(L), Ell}, {x, y},explicit);
min(map(t -> eval(y, t), [%]));
MySol := select(s -> has(s, y=%), [sol])[];
point(N2, eval([x, y], MySol));
segment(sg, N1, N2);
tx := plots:-textplot([[coordinates(M)[], "M"],[coordinates(N1)[], "N1"],[coordinates(N2)[], "N2"],[coordinates(P)[], "P"],[coordinates(Q)[], "Q"],[coordinates(F)[], "F point de Frégier"],[coordinates(T)[], "T"]], font = [times, bold, 16], align = [above, left]);
dr := draw([sg(color = magenta, linestyle = dash),Per(color = black), P(color = red, symbol = solidcircle, symbolsize = 12),Q(color = red, symbol = solidcircle, symbolsize = 12),M(color = black, symbol = solidcircle, symbolsize = 12),F(color = red, symbol = solidcircle, symbolsize = 12),N1(color = black, symbol = solidcircle, symbolsize = 8)]);
plots:-display(Courbe, tx, dr, scaling = constrained, axes=normal, view=[-15..15, default]);
end proc:
plots:-animate('OneFrame', [k], k = Pi/3 .. Pi, frames = 50);
 

(mw file on demand, in particular if you want to display correctly french accents)

Animation 

@acer 

I never intended to present something concise, not even to spek of delivering a definite answer.
My goal was to make @WA573 understand, through a step by step approach, how to manipulate an expression containing RootOf, why there was not a single solution, and why thecontent of ots last replay to you was inconsistent.

You can delete my comment if it doesn't suit you or if you find it so “unnecessarily complicated” that it only confuses the OP.

@C_R 

, this "However, students look for fun and cool stuff," phrase.

And I ine hundred percent agree with you saying "I am reluctant to waste my spare time on CAD tools since I only do these simulation experiments for fun whenever I need distraction from boring duties."
For years I've thought it was a huge waste of time and that only the quality/reliability on the results mattered... until I presented my works to decision makers who, at the image of your students, only look for fun and cool stuff.

So, instead of showing them some displacements and forces versus time curves which were, supposedely, two "abstract", I spent hours to do that (where I oversimplified a lot of things, in particular about the true movement of the gold balls).

@WA573 

but you have to give  acer some clarifications:

Need_clarification.mw

As I titled I don't want to interfere between you and acer,  so let me know once you haved read this comment in order I might delete it.

@Ronan 

Do you mean something like that?

example_2.mw

@Ronan 

Didn't you already asked almost the same question here one and a half month ago?

What @Christopher2222 is suggesting seems close to the idea I developped when I answered your initial answer.

example.mw

@C_R 

I should be the one to thank you for having taken the time to reply to these technical questions.
I'm really blown away by what MapleSim can do now, and those sphere cylinder contact elements are just those I needed  a few ears ago.
I hope I'll see you soon with more creations of this kind.

By the way, in case you don't already know Tadashi Tokieda's works, here is a Youtube video where it presents some physical paradoxes... new challenges for MapleSim (and you)?

@C_R 

I indeed wondered if it could be a visual artefact.

I didn't carefully read all the text of your post (my apologies) and I have a question to which you may have already given some answers. 
You mention frictional effects and I feel that the movement of the ball within the tube is governed by rolling (where no dissipation occurs) and slipping (which dissipates energy).
Is MapleSim capable of determining by itself whether the ball slips and rolls, slips without rolling, rolls without slipping or it is you who chose a priori what is about to occur?

A suggestion: Could we draw an equator and meridians on the ball to see its movement?

I'm asking this question because six or seven years ago, I worked on an industrial mechanism where a rigid ball was squeezed between two rigid surfaces with relative displacement. I did a complete study of this the behaviour of this ball using Hertz's contact theory and, as the final equations were quite complex, I thought of using MapleSim (which I was using at the time) to model and simulate the problem.
But either I wasn't competent enough or MapleSim didn't have the capabilities it should have today, so I eventually gave up.

By the way, you did (MapleSim too :-) ) quite an impressive work!