Very nice.
Two comments on the worksheet itself:
1. You are apparently setting kernelopts(opaquemodules=false). Otherwise `pdsolve/BC`:-EvaluateBCAtSolution wouldn't work.

2. You have a line which in my downloaded version looks like simplify( Equality - Equality), which to the uninitiated obviously is zero and should be without simplify. Apparently some labelling is going on which I don't see, probably because I don't use labels ever.

@mmcdara Actually the result of dsolve(sys) doesn't look bad if the RootOfs are replaced by aliases (or just by names).
It seems that Im can decide if the roots areal.
 

restart;
sys := {
diff(tau__1(t),t)=phi__1(t),
diff(tau__2(t),t)=phi__2(t),
diff(phi__1(t),t)+2*phi__1(t)+3*phi__2(t)+4*tau__1(t)+5*tau__2(t)=6,
diff(phi__2(t),t)+8*phi__1(t)+7*phi__2(t)+9*tau__1(t)+10*tau__2(t)=11
};
sol:=dsolve(sys);
RO:=indets(sol,specfunc(RootOf));
Im~(RO); # {0} so all roots are real
Re~(RO); # Just returns the input consistent with all roots being real

I don't use alias very often, but here it might be useful:
 

alias(seq(alpha[i]=RO[i],i=1..4));
sol; # Short
evalf(sol);

Be aware of the following point in the help for alias:

"Because aliases are resolved at the time of parsing the original input, they substitute literally, without regard to values of variables and expressions.  For example, alias(a[1]=exp(1)) followed by evalf(a[1]) will replace a[1] with exp(1) to give the result 2.718281828, but evalf(a[i]) will not substitute a[i] for its alias even when i=1.  This is because a[i] does not literally match a[1]."

This implies that the following attempt to find floating point values for the roots won't work:
 

seq(alpha[i],i=1..4);
evalf(%); # no floats, just names
## The simple way works
RO;
evalf(%);

Ordinary names could be used instead of aliases as done here:
 

alias(seq(alpha[i]=alpha[i],i=1..4)); #Removing the aliases first
RO; 
#sol; # Now big
sol2:=subs(seq(RO[i]=alpha[i],i=1..4),sol); 

sol2 is as short as sol with alias, but of course evalf won't do anything about the roots in sol2.

A last comment: I tried evalc  in order to get an expression that also on the surface looks real:
 

r1:=allvalues(RootOf(_Z^4+9*_Z^3+4*_Z^2-19*_Z-5, index = 1));
r1a:=evalc(r1);

I succeeded in some attempts (after restarting) and other times it took too long for my patience.
By using indets(r1a,function) I found these:
{arctan(3/3268*sqrt(10275765)), cos((1/3)*arctan(3/3268*sqrt(10275765))), cos(2/3*arctan(3/3268*sqrt(10275765))), sin((1/3)*arctan(3/3268*sqrt(10275765))), sin(2/3*arctan(3/3268*sqrt(10275765)))}
All very good, but indets(r1a,radical) brought some square roots of not oviously positive quantities.

@annarita To get that error you must be using 2D input.
In 2D input a space is interpreted as multiplication. My code had a space between RGB and ( R() , R(),R() ).
While that is fine in 1D input (aka Maple input) it isn't in 2D input.
Just remove the space so that you have
 

 RGB( R() , R(), R() )

I tested your code in Maple 12 and Maple 15 besides in Maple 2018 and didn't run into any problems. I don't have Maple 13 on this machine.
The 3 results agreed after simplification.

@Carl Love odetest works fine on the implicitly given solutions regardless of the names or values of the constant.
On the explicit solutions having the RootOf odetest converts to an implicit form, which is seen by setting infolevel[odetest] to 2 or higher:
 

restart; 
ode := diff(y(x), x) = x*ln(y(x)):
sol:=dsolve(ode, implicit);
odetest(sol,ode); # 0
odetest(subs(_C1=C1,sol),ode); # 0
odetest(subs(_C1=1,sol),ode); # 0
explicit_sol:=op(solve(sol,{y(x)}));
infolevel[odetest]:=2:
odetest(explicit_sol,ode); # 0
odetest(subs(_C1=C1,explicit_sol),ode); # problem
odetest(subs(_C1=1,explicit_sol),ode); # problem

The procedure used for the job is `odetest/implicit`.

@pppc I used Q as a constant and should have mentioned it above.
The method still works for Q = (1/2)*q*(x-50)^2-1250*q, however, but you need to give q a value before plotting, e.g. q=1 as in

plot(eval(U, {C__1 = 0.12e9, C__2 = 0.2e10, C__3 = 0.2e9, C__4 = 0, C__5 = 0.3e8, q = 1}), x = 0 .. 100, size = [2000, default]);

@guto_vasques You simply forgot to write range = ... . Surely you meant to write:
 

p2b := dsolve(sys2, type = numeric, abserr = 1.*10^(-12), relerr = 1.*10^(-12), range=te .. tr);

@Carl Love In the help page dsolve, numeric, DAE we find the statement:
"The following options are also available for some or all of the DAE methods:
....
implicit = boolean
... "
I tried the first example given on the help page. It uses the methods rkf45_dae (the default), rosenbrock_dae, and mebdfi.
For the first two of these the option implicit=true was accepted and gave slightly different results than without thereby indicating that a different computation was made. With method = mebdfi and implicit=true the error message was

Error, (in dsolve/numeric) 'implicit' can only be used with the rkf45, ck45, rosenbrock, dverk78, lsode, classical or gear methods

@sherek The imaginary unit in Maple is by default named I (capital). You are using the letter i instead.
If you try
 

interface(imaginaryunit);

and the answer is I, then you have the default setting.
You can change to another name such as j if you wish (as some people in electical engineering do).
Just do
 

interface(imaginaryunit=j);

The response this time is the previous setting I. Don't let that confuse you.

As always an example of this error happening will help you get an answer.

It is easy enough to produce the error if you feed simplify with some nonsense as in this example:
 

simplify(xxx,a=yyy);

where xxx, yyy, and a are unassigned.
But I may assume you are not talking about nonsense like that?

@Christopher2222 As we agree (I think) the ball stops bouncing in a finite time.
When I said that if we don't stop at some finite number of hits, then dsolve/numeric will do it for us I was referring to that error/warning message. Also there is roundoff/truncation or whatever to take into account.
That the integration stops because of Event2 could be considered a more graceful stop than a stop due to maxfun exceeded, but it is probably not worth the trouble.

If we try 1000 hits and raise maxfun nothing much is obtained. Try e.g. maxfun=10^6. Then Event2 is triggered at 
t = 5.19806251391305  (try dsol(5.2)  ) .
Not a whole lot different from the setting of 100 hits, where Event2 is triggered at
t = 5.19790291862644. (again try dsol(5.2) ).
Note: If we try lowering abserr and relerr we easily run into warnings about a singularity.

@Christopher2222 You have to stop the infinitely many events from happening or dsolve/numeric will do it for you.
After that it becomes another ball game, which is that of a rolling ball constrained to move on a given surface.
This modified version of Rouben's worksheet uses a discrete variable b(t) to count the number of hits and stops after 100 hits.
The surface is just a parabola.
MaplePrimes18-04-21_bouncing-ball-on-hilly-terrain.mw

@Rouben Rostamian  Does the total energy decrease in your revised model?
Shouldn't it?
OK, I see that it appears so:

odeplot(dsol, [t,(diff(x(t),t)^2+diff(y(t),t)^2)/2+y(t)*9.81],0..10);

Your problem has already been pointed out by acer:  e is just a name with no value, unless it has been assigned to.
In all but very old versions of Maple a warning message comes up if you try dsolve/numeric on an initial value problem with unassigned parameters. The warning from

res:=dsolve({D(y)(x)=(-2*e^x)/(2+e^x), y(0)=1/9}, type=numeric);

is
Warning, The use of global variables in numerical ODE problems is deprecated, and will be removed in a future release. Use the 'parameters' argument instead (see ?dsolve,numeric,parameters)

at least from Maple 12 and up (and including the present release Maple 2018).

I see that you have Maple V (release number irrelevant)  I checked in Maple 8 and there was no warning message either.

I tried the code above in Maple 8 and in Maple 2018, and after that made the assignment

e:=evalf(exp(1)); # getting e := 2.718281828

Then tried:
 

res(1.2345);
res(-1.2345);
plots[odeplot](res,[x,y(x)],x=-1..2);

All 3 lines worked in Maple 12 and Maple 2018 (but the method certainly cannot be recommended).
The first line worked in Maple 8 and so presumably also in your version.
The other two lines didn't work.
Don't mess with this: Use exp(x).
##Note: I recall that in old versions (older than Maple 8) the name E (capital) stood for exp(1).
I cannot recall when that stopped.
Try in your Maple:

evalf(E);

if 2.7128 ... comes up you have such a version. If E is the answer, you don't.

@Carl Love I checked the following trivial pdes in Maple 12, 15, 16, 17, and 18:
 

eq1:=diff(u(x,t),x)=0;
eq2:=diff(v(x,t),t)=0;
dsolve({eq1,eq2},{u(x,t),v(x,t)});

Maple 12, 15, and 16 responded with
Error, (in dsolve) not an ODE system, please try pdsolve
whereas the more recent ones acted as pdsolve would.

I prefer keeping things separate. It helps us think what we are actually doing. It is part of clarifying the problem before starting to solve it.

Thus I am not pleased with the command PDEtools:-Solve either. It is described this way:

"The Solve command computes the value of solving_variables that solves a system of equations sys. The system being solved can involve algebraic or differential equations, or both. You can request an exact (default), numeric, or series solution (respectively use the option numeric or series). In this sense, Solve is a unified command that understands when to call solve, fsolve, dsolve, or pdsolve according to your input, also facilitating the analysis of different types of solutions by just adding the keywords series or numeric." [my emphasis]

This makes me think about my response to students, who faced with some problem would automatically try solve.
My response: solve( "What is the meaning of life?" );