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.

@gaurav_rs I just tried in some old versions (besides the new Maple 2018):
It works in Maple 8, Maple 12, and Maple 15.
I didn't check any other. So maybe you are just more impatient than me. It surely took less than 1/2 minute in any of those versions. I didn't time it.

## OK now for the fun of it I timed factor(A11) and factor(A11,complex) in Maple 8:
The first took 15.155s the second 14.452.

@Adam Ledger Try this for fun (?):
 

do 999 end do;

You need to read about the loop structure in Maple.

@mehdi jafari Still using indets:
 

restart;
f2:=-2*(sqrt(6*beta[11]^2-8*beta[11]*sigma[11]+12*beta[12]^2-24*beta[12]*sigma[12]+4*sigma[11]^2+12*sigma[12]^2)*(-sigma[11]^2-3*sigma[12]^2-3/2*(beta[11]^2)+2*beta[11]*sigma[11]+6*beta[12]*sigma[12]-3*beta[12]^2)+E*delta_gamma*omega*(beta[11]^2+6*beta[12]^2-12*beta[12]*sigma[12]+6*sigma[12]^2))*(1/sqrt(6*beta[11]^2-8*beta[11]*sigma[11]+12*beta[12]^2-24*beta[12]*sigma[12]+4*sigma[11]^2+12*sigma[12]^2))*(1/(3*beta[11]^2-4*beta[11]*sigma[11]+6*beta[12]^2-12*beta[12]*sigma[12]+2*sigma[11]^2+6*sigma[12]^2))
indets(f2,radical);
rdc:=op~(1,%);
subs(op(rdc)=phi^2,f2);
simplify(%) assuming phi>0;

@9009134 I never had a dog, but I like the terrier's tenaciousness.
Your last link https://www.mapleprimes.com/view.aspx?sf=247475_Answer/youssef2010.pdf has the same problem as the one I pointed out earlier. In eq. (59) the laplace transform has been applied to eq. (56) and using eq. (60) that means that he must be assuming that n = alpha-1 > 0.

Undoubtably I'm missing something, but I don't see what it is.
To see the problem illustrated in Maple ( but considering the actual problem itself: the integral):
 

restart;
H:=(f,alpha,x,t)->Int((t-s)^(alpha-1)*f(x, s)/GAMMA(alpha), s = 0 .. t);
f:=(x,t)->x^2*t; # Example
H(f,0.95-1,x,t); # alpha-1
value(%) assuming t>0; #A serious problem
H(f,0.95,x,t); # alpha
value(%) assuming t>0; # No problem