@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

@Mariusz Iwaniuk The three authors of the paper from which we get equation (16) and (17) don't attempt any explanation of how they got those results.
They go on to say that they have used that (in their notation):
L{I^alpha f(t) } = 1/s^alpha*L{f(t)} for alpha > 0 (their equation (18) ).
In the pde they have I^(alpha-1), thus it appears that they need alpha -1 > 0.
But they take alpha = 0.95, so I'm lost.

@Mariusz Iwaniuk In the reference https://mapleprimes.com/view.aspx?sf=247318_Answer/c6ra28831f.pdf given by the OP just before equation (16) it says:
"The closed-form solution ... in the Laplace domain can be found as: "
and then follows equation (16) which has r1 and r2 as functions of the laplace variable s.
It appears to me that you invert r1 and r2 numerically (and not the total expression).
In other words you take theta( Zeta, beta) to be the expression obtained by replacing r1 and r2 with the inverse laplace transform.
Is that really the meaning of what the authors write?