@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]);

@shrey183 The many curves shown above are orbits of the solutions in phase space: x(t) - values along the x-axis and diff(x(t),t) -values along the y- axis.
The graphs not shown above are graphs of the solutions, i.e. giving x(t) - values along the y-axis and t - values along the x-axis.
But here are those graphs:

Notice that for a given value of t0 and x0 there are infinitely many solutions passing though (t0,x0), one for each slope i.e. for each value y0 of diff(x(t),t) at (t0,x0). Thus it doesn't make any sense to ask for arrows in this plot.

@Math Pi Euler You can do something like this:
 

plots:-display(Array( 
   [seq(
      [seq(plots:-odeplot(Resfu(d),[x,z(x)-eval(Z[k](x),delta=d)],caption=typeset("Error for order ",k, " and delta = ",d )),k=0..3)],
      d=[0.3,0.2,0.1]  )]
   ));

where the result is an array with 3 rows and 4 columns because you got a list of lists of plots like this:
[ [a1,a2,a3,a4], [b1,b2,b3,b4], [c1,c2,c3,c4] ].
MaplePrimes won't accept pasting the array of plots, but you can try it. It worked for me.

I have tried changing the coefficients in the polynomial x^6 -3*x^2-2*x+11 to other integer values and have found no problem when changes were made. I was still using your expression ee.
So a natural question is: What is so special about this particular polynomial?

Just a note on evala.
In versions of Maple from 2016 and earlier your evala(ee) doesn't do anything whatsoever to ee.
If you do

op(3,eval(evala));

in Maple 2016 you get:
`Copyright (c) 1991 by the University of Waterloo. All rights reserved.`
If you do the same in Maple 2017.3 you get
lock, `Copyright (c) 2016 by Maplesoft. All rights reserved.`
and in Maple 2018:
`Copyright (c) 2016 by Maplesoft. All rights reserved.`

I did this many times:
 

restart;
a:=[seq(RootOf(_Z1^6-3*_Z1^2-2*_Z1+11, index = i), i = 1 .. 6)];
ee := a[1]*a[5]+a[2]*a[6]+a[3]*a[4];
evala(ee);

The first time I lost the kernel connection.
Several times I got the correct value 4.
One time I got the alarming error message:

Error, (in evala/Normal/preproc0) this should not happen

I have chosen interface(typesetting=standard) in Tools/Options/Display so don't get the problem with INTERVAL.
This works fine:

interface(typesetting=standard);
INTERVAL(1..2, 3..4);

but not if typesetting=extended is used.

I changed the order of alias and seq and repeated the following code including the restart 10 times in one worksheet.
Then executed the worksheet using !!! and doing it several times. I got the result 3.139680582-1.737673929*I each time except for the last digits which varied between 2, 3, and 4 for the real part and 27, 28, 29, 30 for the imaginary part.
 

restart;
alias(seq(a[i] = RootOf(_Z1^6-3*_Z1^2-2*_Z1+11, index = i), i = 1 .. 6));
ee := a[1]*a[5]+a[2]*a[6]+a[3]*a[4];
evalf(evala(ee));

I also tried inserting an evalf(ee) before the evalf(evala(ee)). That made it very likely that 4. was the result, but not consistently. The other result was as above.
 

restart;
alias(seq(a[i] = RootOf(_Z1^6-3*_Z1^2-2*_Z1+11, index = i), i = 1 .. 6));
ee := a[1]*a[5]+a[2]*a[6]+a[3]*a[4];
evalf(ee);
evalf(evala(ee));

I use the 64 bit version of Maple 2018.
## Addendum
The following also mostly return 4 but not consistently:
 

restart;
alias(seq(a[i] = RootOf(_Z^6-3*_Z^2-2*_Z+11, index = i), i = 1 .. 6));
ee := a[1]*a[5]+a[2]*a[6]+a[3]*a[4];
evala(a[1]); # has an effect apparently
evala(ee);

You cannot have infinity in your boundary conditions.
You may have a look at the building blocks for a possible series solution:
 

pdsolve(pde,build);

@Math Pi Euler Just do Z[3](0) if x = 0.
I added a comparison to the numerical solution as done in my other answer.