@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?" );

Your ball stops at a finite time after an infinite amount of jumps, doesn't it?
Some quick calculation suggests that it stops at t = 3*sqrt(2), which is approximately 4.242640686 thus it may not be so surprising that you got problems. But you were probably aware of that.
## So dsolve/numeric will be expected to handle an infinity of events!  :-)

This brings up the often lamented fact that while the help page for dsolve,numeric,events mentions quite a lot of possibilities, there are only two examples given on that page.
Those two examples don't illustrate (and cannot be expected to illustrate) all the many features apparently available according to the help page. Thus this presumably great work is basically out of reach to the average user.
Documentation is badly needed!

@tomleslie An exact solution (as in y(x) = sin(x)*exp(x) ) is nice to have if the functions appearing there have well documented properties (as do sin and exp). It is an additional advantage if those functions are also well known to the actual user.
An exact solution as in our case involving AiryAi, AiryBi, and unevaluated integrals involving those functions surely could find some use in some situation. In the present case, however, where the only important thing is to grind out numbers for a comparison with an approximate solution (found by using a wavelet approach) I would strongly favor using dsolve/numeric as I have already stated.
After all, the uniquely given maximal solution to the initial value problem

ode:=diff(Y(t),t,t)-diff(Y(t),t)+t*Y(t)=t^2;
ics:= {Y(0)=1,D(Y)(0)=-1});

is defined on all of the real axis and may be given a fancy name, why not Fellini, and the recipe for computing that exact solution Y(t) = Fellini(t) is simple: Use dsolve/numeric.

Even the help page for AiryAi and AiryBi starts by referring to the basic property of those functions:

"The Airy wave functions AiryAi and AiryBi are linearly independent solutions for w in the equation w'' - z*w = 0."

@tomleslie You say " Maple can solve this ODE exactly, so I have used this as the "exact" solution ".
OK, in some sense yes. But the solution is given in terms of integrals that Maple cannot find exactly:
 

ode:=diff(Y(t),t,t)-diff(Y(t),t)+t*Y(t)=t^2;
sol:=dsolve({ode,Y(0)=1,D(Y)(0)=-1});
value(sol);

Y(t) = -(1/2)*exp((1/2)*t)*AiryAi(1/4-t)*(3*AiryBi(1/4)-2*AiryBi(1, 1/4))/(AiryAi(1/4)*AiryBi(1, 1/4)-AiryAi(1, 1/4)*AiryBi(1/4))+(1/2)*exp((1/2)*t)*AiryBi(1/4-t)*(3*AiryAi(1/4)-2*AiryAi(1, 1/4))/(AiryAi(1/4)*AiryBi(1, 1/4)-AiryAi(1, 1/4)*AiryBi(1/4))+exp((1/2)*t)*Pi*(AiryAi(1/4-t)*(int(AiryBi(1/4-_z1)*_z1^2*exp(-(1/2)*_z1), _z1 = 0 .. t))-AiryBi(1/4-t)*(int(AiryAi(1/4-_z1)*_z1^2*exp(-(1/2)*_z1), _z1 = 0 .. t)))

Thus these integrals will have to be computed numerically. The advantage of the "exact solution" is therefore hard to see.

This is not an answer, but just a comment about the production of the data from the worksheet.
Your worksheet doesn't give the data, but it can be generated like this after your worksheet is run:
 

interface(rtablesize=16); #To see some numbers (default is 10)
resM:=dsolve({ode,Y(0)=1,D(Y)(0)=-1},numeric,abserr=1e-15,relerr=1e-13,output=Array([seq((2*i-1)/32,i=1..16)]));
A:=resM[2,1]:
yap:=y~(A[..,1]);
err:=A[..,2]-yap;
M:=<yap|A[..,2]|err>;

Notice that the author of the paper from which your image is taken uses a fifth order polynomial as representing the exact solution. I'm using the result from dsolve/numeric which represents the exact solution much better. I added a comment about that in your other question:
https://www.mapleprimes.com/questions/224390-The-Code-Doesn-T-Work-Why#comment247856

@student_md I wouldn't be the right one to answer that one.
But my guess is that the best way is to export the numbers as a matrix to a text file using either Export or ExportMatrix.
Then do the latex work outside of Maple.
But feel free to ask a new question in this forum about that general problem.
######################################
Note: It appears that the author of the paper uses a truncated series solution using highest degree 5.
That is not particularly accurate and can hardly be claimed to represent the exact solution.
You are much better off using the numerical solution as representing the exact solution.
You could of course increase the order as in this example:
 

resS:=dsolve({ode,Y(0)=1,D(Y)(0)=-1},Y(t),series,order=13);
yS:=convert(rhs(resS),polynom);

But yS is still an approximation and so is the numerical solution, but I suggest using the latter because of the error control available. With the default setting of Digits (10) you can even use abserr=1e-15, relerr=1e-13 as in this:

res:=dsolve({ode,Y(0)=1,D(Y)(0)=-1},numeric,abserr=1e-15,relerr=1e-13);

@student_md I have cleaned the code up some. Removed t as formal parameter from hi and pn since t is being used globally in other places. One thing I haven't done is to replace the indefinite integrals int(hd[i],t) and int(pn(i,n-1),t) by definite ones. You should be aware of that. Are you sure that the integral returned is the one you want?
 

Reference:
http://www.m-hikari.com/ams/ams-2012/ams-125-128-2012/sunmonuAMS125-128-2012.pdf

restart;
h1 := t-> piecewise(0 <= t and t < 1, 1):
###### We will be using t as the global time parameter throughout. #####
hi:=proc(j,k) local a,b,c,m; 
  m:=2^(j);
  a := k/m; 
  b := (k+1/2)/m; 
  c := (k+1)/m;
  return piecewise(a <= t and t < b, 1, b <= t and t < c, -1)
end proc:
##
J:=3: 
N :=2^J:
hd := Vector(N): 
H := Matrix(N, N): 
T := Vector(N):
hd[1] := h1(t):
for i from 1 to N do T[i] := (i-1/2)/N end do:
##
for j from 0 to J-1 do
  m := 2^j;
  for k from 0 to m-1 do
    i := m+k+1;
    hd[i] := hi(j, k)
  end do
end do:
#Now Compute H at the collocation points
for i from 1 to N do
  for j from 1 to N do
    H[i, j] := eval(hd[i], t = T[j])
  end do
end do:
##
pn:=proc(i,n) 
  if n=1 then 
    return int(hd[i],t) ## Notice the use of an indefinite integral!
  else
    return int(pn(i,n-1),t) ## Notice the use of an indefinite integral!
  end if
end proc:
##
##EXAMPLE 2 from the reference given:
f:=t->t^2;
alpha1:=t->-1:
alpha2:=t->t:
y0:=1:
y1:=-1:
##
RHS:= t->f(t)-alpha1(t)*y1-alpha2(t)*(t*y1+y0);
R:=Vector(N):
TMP:=Matrix(N,N):
A:=Matrix(N,N):
##
for i from 1 to N do
  R[i] := evalf(RHS(T[i]));
  tmp := alpha1(t)*pn(i,1)+alpha2(t)*pn(i,2);
  for j from 1 to N do
    TMP[i,j]:=eval(tmp, t = T[j])
  end do
end do:

##Compute the wavelet coefficients:
use LinearAlgebra in
  A := Transpose(LinearSolve(Transpose(H+TMP), R))
end use:
#Now compute the approximate solution
sol := y1*t+y0+add(A[m0]*pn(m0,2), m0 = 1 .. N):
#Convert to a function of t for easy comparison with exact solution
y:=unapply(sol,t):
#############################################
## Compare with the solution found by dsolve/numeric:
## EXAMPLE 2 ode:
ode:=diff(Y(t),t,t)-diff(Y(t),t)+t*Y(t)=t^2; 
res:=dsolve({ode,Y(0)=1,D(Y)(0)=-1},numeric,abserr=1e-10,relerr=1e-9);
plots:-odeplot(res,[t,Y(t)-y(t)],0..1,caption="The error on the interval 0..1");

########## In the text in the reference I see that integrals from 0 to t are used.
In recent versions of Maple (certainly from 12 and up, but not in Maple 8) you can use the physicists' abuse of notation t=0..t: Example:

u:=exp(t):
int(u,t); # not zero at t=0.
int(u,t=0..t); # t is used in two roles.
## Alternatively, without abuse of notation and works in all versions:
subs(tt=t,int(u,t=0..tt));

The best "splash screen" I remember was in Maple V, Release 3. It had a portrait of Isaac Newton.

The image you are showing makes me think about Grønlangkål og stuvet hvidkål, a Danish cabbage dish served mainly at Christmas.
The blue color is distracting though, since that dish is usually served with medisterpølse, which is a sausage and not blue at all. But maybe artistic license should be allowed for.

@rlewis Markiyan Hirnyk was objecting to Kitonum's answer, because you in your question wrote that the real example is much more complex.
Kitonum made f, g, and h into procedures using the arrow notation:
f:=(x,t)->t^2*x^2 + t*x + 2*x - 1;
Notice that if you do that in Maple and copy the blue output and paste it into an input line you get
f := proc (x, t) options operator, arrow; t^2*x^2+t*x+2*x-1 end proc
The first way of writing f is very convenient and is understood by the parser. But it is exactly the same procedure as the second version. In fact if you press enter after putting the second version on an input line then it prints like the first version because of option operator, arrow.

@rlewis No, I mean this:
 

ics:= { x(0) = x0, y(0) = y0, z(0) = z0};

for the correct values of x0, y0, and z0.

@rlewis OK, if you already know the initial point then just use that in dsolve/numeric.
In other words, just skip the lines

eval(eqs2,t=0);
ics:=solve(%,{x,y,z}(0));

Then define odes as above and write down the correct initial point in the same form as ics above.