@mmcdara Notice that before dsolve sees the ode the right side with `if`:
`if`(t=1, 0, x(t));
evaluates to x(t) for the trivial reason that the name t is not equal to 1.
That is why you need piecewise.

@radaar I don't quite understand what you mean.

But incidentally, I found a bug in dsolve/symbolic:
 

restart;
ode:=diff(x(t),t)=piecewise(t=1,0,x(t));
sol:=dsolve({ode,x(0)=1});
odetest(sol,ode);
plot(rhs(sol),t=0..2,caption="The symbolic solution is wrong");
res:=dsolve({ode,x(0)=1},numeric);
plots:-odeplot(res,[t,x(t)],0..2,caption="The numerical solution is correct");

The symbolic solution is wrong, the numerical one is correct.

@torabi Well, you didn't really answer my question about the amplitude: What does that mean when the attractor is rather complicated and not just a point?
So what I decided to do was to sample the y-values for various time values equidistantly spaced after an initial settling period.
The code follows.
 

## This assumes that you have executed the definition of my version of fdsolve given earlier.
## option remember in the procedure means that if you repeat any of this for the same
## A-value then it won't take long. 
##That is why I stick to a previously defined list of A-values.
##
F:=omega*x-y^2;
G:=mu*(z-y);
H:=A*y-B*z+x*y;
params:={omega=-2.667,mu=10,B=1}; # A unassigned
## Order alpha, final time T, number of points N. All will be fixed in advance.
alpha:=0.89: T:=20: N:=2000:
solfu:=proc(AA,inits::[numeric,numeric,numeric]) option remember,system; local FGH;
    if not AA::numeric then return 'procname(_passed)' end if;
    FGH:=unapply~(eval([F,G,H],params union {A=AA}),t,x,y,z);
    fdsolve(FGH,alpha,0..T,inits,N)
end proc:
## We shall be using two initial conditions inits1 and inits2:
inits1:=[.1,.1,.1]; inits2:=[-.1,-.1,-.1];
## The A values considered are:
Alist:=[seq(0..30,1)];
sol:=solfu(Alist[1],inits1); #Test (not wasted because of option remember)
plot(sol[..,[1,3]],labels=[t,y],size=[1200,default]);
## The following animation is not wasted because of option remember.
plots:-animate(plots:-spacecurve,[solfu(A,inits1)[..,2..4] ],A=Alist);
## For each A we pick points beginning with Nb and continuing with a spacing of n points.
## Thus we pick the points numbered Nb+i*n, i=0..floor((N-Nb)/n).
Nb:=N-1000; n:=50;
numpts:=floor((N-Nb)/n)+1; # Number of points chosen for each A 
## The distance in time between points will be
n*T/N;
## The procedure Q produces a sequence of points for A given.
Q:=proc(A,inits::[numeric,numeric,numeric]) local sol;
   sol:=solfu(A,inits);
   seq([A,sol[Nb+i*n,3]],i=0..numpts-1)
end proc;
##
PtList1:=[seq(Q(A,inits1),A=Alist)]: # This is fast because of the animation done earlier.
plot(PtList1,style=point); p1:=%:
PtList2:=[seq(Q(A,inits2),A=Alist)]: # This takes a while.
plot(PtList2,style=point,color=blue); p2:=%:
plots:-display(p1,p2,labels=[A,y],size=[800,700],caption=typeset('alpha'=alpha));

The plot:

 

Because of option remember the following won't take long:
 

## It shows animations in A of the orbits from inits1 and inits2 in the same coordinate system.
use plots in
  a1:=animate(spacecurve,[solfu(A,inits1)[..,2..4] ],A=Alist);
  a2:=animate(spacecurve,[solfu(A,inits2)[..,2..4] ],A=Alist);
  display(a1,a2)
end use;

@Mariusz Iwaniuk A detail about Pi which you have outside evalf and which would cause a problem (easily cured though) in earlier versions of Maple.
In versions from Maple 18 and earlier the option kernelopts(floatPi) didn't exist. The default setting in the Maple releases 2018, 2017, and 2016 is floatPi=true. The option is available in Maple 2015, but is not documented; default is floatPi=false.
What it means having kernelopts(floatPi=true)  is that the presence of a float in a product or sum involving Pi makes Pi evaluate as a float.
 

kernelopts(floatPi=true);
Pi*1.0; #Results in a float
Pi+0.0;  #Results in a float
[Pi,1.0]; # Doesn't

Personally I don't like that setting and have kernelopts(floatPi=false) in my maple.ini file.
One reason for disliking the setting is that it singles out the constant Pi. It doesn't apply to any other constant.
Example: sqrt(2)*1.0 evaluates to just that, not to a float.

@Mariusz Iwaniuk The error is not impressively low:
 

plot(sinh(x)-ifunc(5)(x), x = -3 .. 3);

It doesn't appear to improve by using more iterations.
I think you would be better off using finite elements with collocation, i.e. replacing y by a linear combination of basis functions (e.g. monomials) Y with coefficients to be determined. Then evaluate the resulting equation at some points and solve for the unknown coefficients.
In this linear case that approach will be far better.

@torabi Plotting y as a function of A assumes that you have a precise definition of what y is.
You stated earlier that y meant " the amplitude of the y component of the system's attractor".
What does that mean when the attractor is rather complicated and not just a point?

@torabi No I haven't looked much at it.
What I have done is to look at the dependence on A for fixed alpha = 0.89.

Using the version of fdsolve given in the worksheet FracDiff.mw and using the example in the paper, I did the following.
 

## Assuming fdsolve defined as in FracDiff.mw
F:=omega*x-y^2;
G:=mu*(z-y);
H:=A*y-B*z+x*y;
paramsA:={omega=-2.667,mu=10,A=AA,B=1}; # A and AA unassigned
params2:=unapply(paramsA, AA); 
## Using alpha=0.89:
##
solfu2:=proc(AA) option remember,system; local FGH;
    if not AA::numeric then return 'procname(_passed)' end if;
    FGH:=unapply~(eval([F,G,H],params2(AA)),t,x,y,z);
    fdsolve(FGH,0.89,0..20,[.1,.1,.1],2000)
end proc;
##
sol:=solfu2(27.3); #Test
plot(sol[..,[1,3]],labels=[t,y],size=[1200,default]);
plot(sol[..,[1,3]],y=10..20,labels=[t,y],size=[1200,default]);
plots:-animate(plots:-spacecurve,[solfu2(A)[..,2..4] ],A=[15,20,25,27.3]);
## The next will be fast because solfu2 has option remember:
plots:-animate(plot,[solfu2(A)[..,[1,3]] ,t=10..20],A=[15,20,25,27.3],size=[1200,default]);

solfu2(A) just computes the solution for a given value of A.
The plots found are just the orbits of graphs of the solutions.
When t starts at a late value (e.g. as here 10) then perhaps enough time has passed for the orbit to be close to a possible attractor.

@Seb1123 Stating the obvious (even without having to consult the help):  simplify will do its best to simplify.
Simplification of an equation just means simplifying each side.
In fact you don't need simplify to do that job for your equation. Just write it down and press ENTER:
F*R=(1/2*M*R^2)*(a/R);  
The output is F*R = (1/2)*M*R*a.
So you expected simplify to read your mind. To read that you wanted to solve the equation for F?

@Federiko My question was basically meant as 'rhetorical'.

If a locally Lebesgue integrable function f defined on [0, infinity[  satisfying for some M >0, some real a, and some t0 >= 0 the inequality
abs(f(t)) <= M*exp(a*t) for a.e. t >=t0  ( "f is exponentially bounded")

then the Laplace transform F(s) exists for all s satisfying Re(s) > a. Furthermore F is analytic in that region. In particular continuous, thus cannot have singularities.
If your F has infinitely many singularities (s[n]) and s[n] -> infinity, then your F cannot be the Laplace transform of an exponentially bounded locally Lebesgue integrable function.

@Federiko In earlier versions of Maple you must replace the type specfunc(sin) with specfunc(anything,sin).
So the code is this (and that will also work in more recent versions):
 

restart;
F:=1/(s^2*(1+(1-1/(2+2*s))/s)*(1+tan((1/2)*Pi/sqrt((1-1/(2+2*s))/s))/sqrt((1-1/(2+2*s))/s)));
F2:=normal(convert(F,sincos));
A:=op([1,1],indets(F2,specfunc(anything,sin)))/Pi; ## Change made here
limit(A,s=infinity);
sp1,sp2:=solve(A=p,s);
evalf(seq(sp1,p=1..10));
dF2:=denom(F2);
seq(fsolve(dF2=0, s=eval(sp1,p=n)..eval(sp1,p=n+1)),n=0..9);

The 'anything' in specfunc(anything,sin) means that the argument to sin must be of type 'anything'. As the name 'anything' indicates this means basically anything. See ?type,anything

@ The code works as written in Maple 2018.1. I just checked the downloaded version from MaplePrimes using the icon !!!
I copied the code and tried it in Maple 18.02.
You are right. In that version parameters appearing in the initial conditions apparently must be "alone", i.e. the parameter H cannot appear as in Jir( H- He ) = ...
## To see a simple example of that, try
 

ode:=diff(x(t),t)=x(t);
res:=dsolve({ode,x(t0-2)=1},numeric,parameters=[t0]); # Error in Maple 18
###
ode:=diff(x(t),t)=x(t);
res:=dsolve({ode,x(t0)=1},numeric,parameters=[t0]); # t0 appearing alone: OK also in Maple 18
res(parameters=[0]);
res(1);

Both versions work in Maple 2018.1:
The solution for Maple 18 is simple. Replace the parameter H by the parameter HmHe (meaning H minus He).
In the initial conditions HmHe shall replace H-He. When you set the parameter HmHe using some value for H, HH say, then you do
resNew( parameters=[HH-He]);
I have uploaded a version with these changes. I have kept the old lines for comparison, but they are commented out.

MaplePrimes18-08-21_odesys_inits18.mw

I changed your code just slightly and didn't use optimize.
The change was that sol was made a function of H defined in advance.
Your version with nn:=20000 executed in 65.5s and mine in 43.7s.
If I add optimize to my version it takes the same as yours.
## The following note is edited: Since the loop breaks at i = 742 you are optimizing dsolve not 20000 times, but only 742 times. That appears to take 20s.
If you remove optimize from your version it takes the same time as mine.
So better ideas are needed.

Now the code anyway:
 

## Before the loop and for clarity unassigning H first (it was a function of t) put this:
##
sys_ode:=diff(sv(x),x)=rho0*g/Lambda(x),diff(sh(x),x)=beta12(x)/beta22(x)*rho0*g/Lambda(x),diff(Lambda(x),x)=rho0*g/beta22(x),diff(u(x),x)=-Mp*g/beta22(x),diff(Jir(x),x)/Jir(x)=-Gp(x)*diff(Lambda(x),x)/Lambda(x);
ics:=sv(H-He)=sve,sh(H-He)=she,Lambda(H-He)=Le,u(0)=0,Jir(H-He)=Jire;
## Now the function of H:
res:=HH->dsolve(eval({sys_ode,ics},H=HH),numeric,[sv(x),sh(x),Lambda(x),u(x),Jir(x)],output=listprocedure);
## Now the loop
nn:=20000:
t0:=time():
dt:=(Ts-Te)/nn:
t:=Te:
HH:=He: #Note HH
for i from 1 to nn+1 by 1 do
 sol:=res(HH);
 CC2(i,1):=t/(10^6*365*24*60*60):
 CC2(i,2):=HH: # HH
 sh0:=rhs(sol(0)[3]):
 sv0:=rhs(sol(0)[2]):
 TR2(i,1):=-1/3*sh0-2/3*sv0:
 TR2(i,2):=sqrt(3)/3*(sh0-sv0):
 Jir0:=rhs(sol(0)[6]):
 phix0:=1-(1-phi0)/Jir0:
 EN2(i,1):=pc0*(ln(phix0)/ln(phi0))^m:
 up:=rhs(sol(HH-He)[5]): #HH
 uh:=ue+up:
 fvp:=evalf(A*sh0+B*sv0-(pv0*(ln(phix0)/ln(phi0))^n)):
 if fvp>=0 then
  break
 else
  dH:=dt*(Mp/rho0+uh):
  HH:=HH+dH: ## HH
  t:=t+dt:
 end if:
end do:
time()-t0;

Maple cannot solve integral equations (or integro-differential equations) numerically. So you are on your own.

First of all you cannot use { }  (or [] ) instead of parentheses ().
Secondly, Maple can solve delay odes, but not delay pdes. In your second equation you have w evaluated at (t - tau, x), where tau = 0.3.

@ecterrab Yes Edgardo, I participated in the discussion in MaplePrimes about matrix input a couple of years ago and thought it would be a good idea. I still think it is. I was primarily thinking of its use in the classroom, where e.g. a system of first order odes with constant coefficients theoretically are handled best as a vector equation x'(t) = A.x(t) + b(t) and handled like that in textbooks.
So thank you much for that feature and for your detailed reply.