@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?
 

@9009134 The laplace transform of the integral appearing on the right hand side of the pde can be found by Maple.
 

inttrans[laplace](Int(theta(xi, s)*(beta-s)^(alpha-2), s = 0 .. beta),beta,z) assuming alpha>1;

Your alpha is 0.95 and as I have already mentioned that is a problem
Even if you consider an alpha > 1 you have to deal with the problem of finding the laplace transform of a product of the integral and another function (one of them being diff(theta(xi, beta), xi, xi) ).
Even if you succeed in getting an expression for the laplace transform of theta(xi,beta) w.r.t. beta from the pde (as the authors of the paper in your link claims to have in (16) and (17) ) then you must invert that laplace transform numerically or by some other approximation method.
Good luck!

@figgisfiddis I just ran your worksheet in Maple 2018 and also in Maple 18.02. I encountered no problem.
I got 0 for the integral in both versions. I use 64 bit Maple on Windows 10.

@9009134 It is a good idea to look at a version with no integrals as you do. It will show you that there are some problems.
First of all notice that there is a syntax difference between the exact pdsolve and the numeric pdsolve.
In the exact mode of pdsolve the pde with conditions are entered together in a set as the first argument.
In the numeric mode the first argument consists of the pde(s) and the second should be a set with all conditions.
Next notice that your pde is first order in xi and 4th order in beta (=time). Thus you can only give one boundary condition, but must give 4 initial conditions.
So I omit one boundary condition here:
 

pdsolve(PDE, Init union Bdry[1..1], numeric,time=beta,range=0..10);

That gives us the error expected from what I said above:

Error, (in pdsolve/numeric/par_hyp) Incorrect number of initial conditions, expected 4, got 2
#############################################################################
The integral

Int(theta(xi, s)/(beta-s)^1.05, s = 0 .. beta);

which appears three times in the original version is divergent if theta(xi, beta) <> 0, which is just because 1.05 >= 1 as is basically shown here:
 

B:=Int(1/(beta-s)^1.05, s = beta0 .. beta);
value(B) assuming beta0>0,beta>beta0;

Result Float(infinity)  (the "Float" is due to the appearance of the floating point number 1.05).
I don't see how using the laplace transform gan get us out that problem.

@Rouben Rostamian  I'm just guessing, but instead of f(t[n+1], X[l-1], X[l-1]) in the code for fdsolve you probably want

f(t[n+1], X[l-1], Y[l-1])

and similarly for g.

@torabi With the change in G and using Rouben's worksheet, but first not his procedure fdsolve, I tried the ode system (not fractional):
 

solve(eval({F(t,x,y)=0,G(t,x,y)=0},params),{x,y});
pts:=map(subs,{%},[x,y]);
J:=unapply(VectorCalculus:-Jacobian(eval([F(t,x,y),G(t,x,y)],params),[x,y]),x,y):
LinearAlgebra:-Eigenvalues(J(op(pts[3])));
sys:={diff(x(t),t)=F(t,x(t),y(t)),diff(y(t),t)=G(t,x(t),y(t))};
ics:={x(0)=0.5, y(0)=0.05};
res:=dsolve(eval(sys,params) union ics,numeric);
plots:-odeplot(res,[x(t),y(t)],4000..4130);

To begin to see the appearance of a limit cycle I went to t=4000..4130.
I started close to the repelling equilibrium {x = .5850885314, y = 0.6038345764e-1}.


After that I tried Rouben's fdsolve with the correction mentioned in my reply to Rouben:
 

T := 1000.0;  # solve over t=0 to t=T
# solution starting at (0.6, 0.35) 
sol := fdsolve(alpha=0.98, t0=0.0, t1=T, x0=0.6, y0=0.35, N=1000, params);
plot([seq([sol[i][2], sol[i][3]], i=0..1000)]);

I got this plot:

@rahinui Since g is unknown (presumably just evaluates to its own name) then why not just do:
 

restart;
f:=(t1,t2)->cos(g(t1,t2));
D[2](f);
## or
D[2](f)(t1,t2);

After all, unapply is used when the direct method used here is impossible because evaluation is needed as in my example given above with gx. The formal parameters t1 and t2 are in plain sight in cos(g(t1,t2)), so no need for unapply.
###
I do find this behavior unfortunate (at least):
 

restart;
f:=cos@g;
D[2](f);

Error, (in D[2]) index out of range: cos takes only 1 argument(s)

Workaround in your case (basically a joke, but it works):
 

restart;
f:=cos@g;
D[1](f)(t1,t2);
eval(%,D[1]=D[2]);