It is interesting that this simple change makes is work as expected
 

is(subs(Pi*x/T=y,eq));

Even this version works:
 

is(subs(Pi=nm,eq));

and therefore also this (just for fun):
 

is(subs(Pi=pi,eq));

The problem in this case seems to be the constant (here Pi). sqrt(2) does no better, but 22/7 works fine.
Even replacing T with 1 gives the wrong result:
 

is(subs(T=1,eq));

###### From a cursory look at the code for `is/solve`  it seems that the problem may be related to the output from testeq, but notice that testeq doesn't return false, but FAIL for eq:
 

restart;
eq := sin(Pi*x/T)*cos(Pi*x/T) = 1/2*(sin(Pi*x/T+Pi*x/T)+sin(Pi*x/T-Pi*x/T));
testeq(eq); # FAIL
testeq(subs(Pi=pi,eq)); # true

A final note: Don't use is when you don't have to, simplifying lhs-rhs is much better here:
 

simplify((lhs-rhs)(eq));

I forgot that you are using Maple 12. There it appears that you need to use expand instead simplify.

1.The fractional powers in the cbs could cause problems, so I would put abs around f(0) where f(0) is raised to the powers 1/3 and 4/3.
 Another fractional power is present in (f(0)^2+(1-g(0))^2)^(1/3), but that should be fine since f(0)^2+(1-g(0))^2 ought to be nonnegative.

2. Moreover I would replace the one occurrence of a dot (.) in bcs with multiplication (*) as is clearly intended.

3. There doesn't seem to be any reason to use the continuation in d you are using. So I would just set d = 1 right away.

4. Finally the approximate solution generated by dsolve/numeric/bvp is too simple. It appears to be:
    [f(eta) = 0., g(eta) = 0, h(eta) = 0, p(eta) = 0, theta(eta) = -eta+1]
I have replaced that with a less trivial one see in the code below.
I first handle one case (beta=0.05,lambda=0.5), then procede with the loops.
 

restart;
with(plots):
Pr := 0.1e-1;
Eq1 := (101-100*d)*(diff(h(eta), eta))+2*f(eta) = 0; Eq2 := (101-100*d)*(diff(f(eta), eta$2))-h(eta)*(diff(f(eta), eta))-f(eta)^2+g(eta)^2-beta*(f(eta)+(1/2)*eta*(diff(f(eta), eta))) = 0; Eq3 := diff(g(eta), eta$2)-h(eta)*(diff(g(eta), eta))-2*f(eta)*g(eta)-beta*(g(eta)+(1/2)*eta*(diff(g(eta), eta))) = 0; Eq4 := diff(p(eta), eta)-2*(diff(f(eta), eta))+beta*(h(eta)+eta*(diff(h(eta), eta))) = 0; Eq5 := diff(theta(eta), eta$2)-Pr*h(eta)*(diff(theta(eta), eta))-Pr*beta*(3*theta(eta)+eta*(diff(theta(eta), eta))) = 0;
##
Vlambda:= [.5, 1, 1.5]; Vbeta:= [0.5e-1, .1, .2, .4, .7, 1, 1.50, 2];
etainf := 1;
d:=1:  # No need for continuation in d 
## NOTICE abs:
bcs:= h(0) = 0, p(0) = 0, theta(0) = 1, D(f)(0) = lambda*abs(f(0))^(4/3)/(f(0)^2+(1-g(0))^2)^(1/3), D(g)(0) = -lambda*abs(f(0))^(1/3)*(1-g(0))* (1/(f(0)^2+(1-g(0))^2)^(1/3)), f(etainf) = 0, g(etainf) = 0, theta(etainf) = 0;
##
dsys:={Eq1,Eq2,Eq3,Eq4,Eq5,bcs};
## FIRST CASE:
IP1:=[f(eta) = 1-eta^2, g(eta) = eta^2/2, h(eta) = 0, p(eta) = 0, theta(eta) = 1-eta];
res := dsolve(eval(dsys,{beta=0.05,lambda=0.5}), numeric,approxsoln=IP1);
odeplot(res,[eta,f(eta)]);
odeplot(res,[eta,g(eta)]);
odeplot(res,[eta,theta(eta)]);
odeplot(res,[eta,h(eta)]);
odeplot(res,[eta,p(eta)]);
############## The loops:
IP:=IP1;
for j from 1 to 3 do   
  for i from 1 to 8 do 
    dsol[j][i]:=dsolve(eval(dsys,{lambda=Vlambda[j],beta=Vbeta[i]}),numeric,approxsoln=IP);
    IP:=dsol[j][i]; # Used in the next i-iteration
    print(i,j); #Just to see progress
  end do;
  ## For the j iteration we revert to the old IP:
  IP:=IP1;
end do;

You could compute fracdiff(x^2,x,t) first as in:
 

res1:=fracdiff(x^2, x, t) assuming t<1,t>0;
res2:=fracdiff(x^2, x, t) assuming t<2,t>1;
plot3d(res1,x=0..1,t=0..2);

Notice that res1 and res2 agree.

You should be differentiating w.r.t. t, not x.

de := 100*diff(v(t),t)=100*9.81-k*v(t);                 
res:=dsolve({de,v(0)=0});
V:=rhs(res);
Vfinal:=limit(V,t=infinity) assuming k>0;
k0:=solve(Vfinal=5,k);
plot([eval(V,k=k0),5],t=0..4);

Just make j start from 2:
 

for j from 2 to N do

I can confirm that the answer to
solve({1<x or x<3 or x<5});
is not satisfactory in Maple 16.02. I get what you get.
In Maple 15 it is the same. In Maple 17, 18, and Maple 2017 two answers come up, one of which is {x=x}. The other {x<5} is not wrong because of course {1<x or x<3 or x<5} is satisfied if x<5.
You could try this
 

is(1<x or x<3 or x<5) assuming x>=5;

The answer is true in Maple 16.

As John Fredsted shows, the ode can be solved exactly (i.e. you get a formula for the solution), but apparently you are asked to solve the ode numerically. So here is a way to do that:
 

restart;
params:={M = 1,G = 1,h = 1};
ode  := diff(1/r(phi),phi$2) + 1/r(phi) - G*M/h^2;
res := dsolve({eval(ode,params),r(0)= 2/3,D(r)(0) = 0},numeric);
## Plotting
plots:-odeplot(res,[r(phi)*cos(phi),r(phi)*sin(phi)],0..2*Pi,scaling=constrained,labels=[x,y]);

I think that Student:-NumericalAnalysis:-InitialValueProblem only works for a single ode.
Use dsolve/numeric instead.
Maybe this is what you want:
 

C:=dsolve({seq(eq[i], i = 1 .. 5), seq(ic[i], i = 1 .. 5)},numeric,method=classical[rk4],stepsize=0.1);

Use Int instead of int and do:
 

A:=(Int((1+t/T)*exp(-I*(2*Pi*n*t/T)), t = -T .. 0)+Int((1-t/T)*exp(-I*(2*Pi*n*t/T)), t = 0 .. T))/T;
IntegrationTools:-Change(A,x=t/T,x);
value(%);
res:=simplify(%);
simplify(res) assuming n::posint;

The last line is only interesting if n is indeed supposed to be a positive integer.

f is a local variable (but z is not). You can use this:

deq2:=convert(deq,`global`);

But notice that since dsolve works without giving the dependent variable, you don't have to do anything but

dsolve(deq);

It would be a problem though if initial or boundary conditions are included.

Your system is of order 3 in F. Thus your boundary conditions can include F, F', and F'', but not F'''.
Your boundary conditions include this one

D(F)(0) = varkappa+delta*(D@@2)(F)(0)+alpha*(D@@3)(F)(0)

To see that that can be a problem try a simple example:

restart;
ode:=diff(y(x),x,x)+y(x)=x;
sol:=dsolve({ode,y(0)=1,(D@@2)(y)(1)=0}); #Exact solution (and correct)
## Now try numerically:
dsolve({ode,y(0)=1,(D@@2)(y)(1)=0},numeric); #Error: range seems to be needed
dsolve({ode,y(0)=1,(D@@2)(y)(1)=0},numeric,range=0..1); # Error: singularity encountered (!)
## We can rewrite (D@@2)(y)(1)=0 by using the ode:
eval(convert(ode,D),x=1);
eval(%,(D@@2)(y)(1)=0); # Thus we have y(1) = 1:
res:=dsolve({ode,y(0)=1,y(1)=1},numeric);
plots:-odeplot(res,[x,y(x)-rhs(sol)]); # Agrees well with the exact solution

You may be able to do something similar in your case.
I strongly suggest making dsolve work for one value in HA before doing a loop.
 

I think the bad news is that there will necessarily be a singularity in the interior when using conditions at x = Pi.
Thus since int(r(y),y=0..x) appears in the original formulation the conclusion is that the problem has no solution.
In the following I use the revised values you gave in your second worksheet.
First the code:
 

restart;
Digits:=15:
with(plots);
R0 := 10^(-5);
C := 1;
m := 1/100;
u1 := 1000;
sys := R(x) = 2*(-r(x)*diff(r(x), x)*cos(x)+r(x)^2*sin(x)+diff(r(x), x)^2*sin(x)-sin(x)*r(x)*diff(r(x),x,x))/(r(x)^4*sin(x)), 
       diff(R(x),x,x)+cot(x)*diff(R(x), x) = -(1/2)*r(x)^2*(R0^2-R(x)^2)-r(x)^2*C^2*m^2*exp(-2*m*A(x))/(2*u1), 
       diff(A(x), x) = r(x);
## I prefer this version:
SYS:=solve({sys},{diff(r(x),x,x),diff(R(x),x,x),diff(A(x),x)});
## I set the condition on A(x) at x = Pi, but of course we don't know its value.
## We let it be a parameter here called APi.
## 
condA := R(Pi) = 1/500^2, r(Pi) = 500, D(R)(Pi) = 0, D(r)(Pi) = 0, A(Pi) = APi;
## Now use the parameters option in dsolve/numeric:
res:=dsolve(SYS union {condA},numeric,parameters=[APi],output=listprocedure);
## Extracting the procedure for finding A(x) for given x:
Aproc:=subs(res,A(x));
## Setting as an example the parameter APi to 10^4:
res(parameters=[10^4]);
## The plots:
plots:-odeplot(res,[x,r(x)],0..Pi);
plots:-odeplot(res,[x,R(x)],0..Pi);
plots:-odeplot(res,[x,A(x)],0..Pi);
## To explore further we use this procedure, which takes APi and x values as input.
## If out=Aproc(xx) is defined for the given APi=aa then that  is the output, otherwise the output is -10^3.
## The most important is what is printed as messages:
##
Q:=proc(aa,xx) local out;
   if not type([aa,xx],list(realcons)) then return 'procname(_passed)' end if;
   res(parameters=[aa]); 
   try 
     out:=Aproc(xx);
     printf("Ok at APi = %f and x = %f\n",aa,xx);
   catch:
     printf(cat(StringTools:-FormatMessage( lastexception[ 2 .. -1 ]),"  APi=%f and x = %f\n"),aa,xx);
     out:=-10^3
   end try;
   out
end proc;
##
Q(2000,1); #Just a single call to Q with APi = 2000 and x = 1.
## Plotting Q(2000,x) (i.e. APi = 2000 for all the x values):
plot(Q(2000,x),x=0..Pi,adaptive=false,numpoints=50);
##
## Now plotting Q(a,1) as a function of a:
plot(Q(a,1),a=1..10^4,adaptive=false,numpoints=50);
plot(Q(a,0.9),a=1..10^4,adaptive=false,numpoints=50);
## It seems to be clear that for APi > 9000 we have reached the maximal interval of definition:
## That interval appears to be 0.99622186 .. Pi.
## It doesn't appear to help to raise APi above 9000.
## The reason is that the exp-term containing A(x) contributes almost nothing for larger values of APi.
##Compare to what happens if we remove that exp-term and therefore also the ode for A(x) as we don't need it:
#Removing the A-part totally from SYS:
SYS0:=remove(has,eval(SYS,exp=0),A(x));
res0:=dsolve(SYS0 union {R(Pi) = 1/500^2, r(Pi) = 500, D(R)(Pi) = 0, D(r)(Pi) = 0},numeric);
odeplot(res0,[x,r(x)],0..Pi); #Notice the value 0.99622186
odeplot(res0,[x,R(x)],0..Pi);

Since you have the integral Ar(x) = int(r(x1), x1 = 0 .. x) you can introduce an ode for Ar:
diff(Ar(x),x) = r(x) with Ar(0) = 0.
You may have a serious singularity at least at r = 0 as will be exhibited in the code below:

restart;
R0 := 10^(-5): C := 1: m := 1: u1 := 1:
sys := R(x) = 2*(-r(x)*diff(r(x), x)*cos(x)+r(x)^2*sin(x)+(diff(r(x), x))^2*sin(x)-sin(x)*r(x)*diff(r(x), x$2))/(r(x)^4*sin(x)), 
      diff(R(x), x$2)+cot(x)*diff(R(x), x) = -(1/2)*r(x)^2*(R0^2-R(x)^2)-r(x)^2*T(x)/(2*u1);
T := x->C^2*m^2*exp(-2*m*Ar(x));
odeAr:=diff(Ar(x),x)=r(x); # Adding this ode
condAr:=Ar(0)=0; # Boundary value for Ar.
cond := R(Pi) = 1/500^2, r(Pi) = 500, D(R)(Pi) = 0, D(r)(Pi) = 0; #Your given boundary values.
##
#Solving the system sys, odeAr for the highest derivatives to see the singularity problem:
SYS:=solve({sys,odeAr},{diff(r(x),x,x),diff(R(x),x,x),diff(Ar(x),x)});
expand(SYS); #Possible problem at x=0 (and perhaps x=Pi)
F:=dsolve(SYS union {cond,condAr},numeric,method=bvp[midrich],initmesh=8192);

I get the error message:

Error, (in dsolve/numeric/bvp) powering may produce overflow

The ode for r(x) in the isolated version is
 

diff(r(x), x, x) = -(1/2)*r(x)^3*R(x)+r(x)-diff(r(x), x)*cos(x)/sin(x)+diff(r(x), x)^2/r(x);

You notice that the last two terms have denominators sin(x) and r(x), respectively. But you need to look at the numerators too, of course.
You could try to avoid initial problems by using an approximate solution where the initial guess for r(x) has r(0)<>0, D(r)(0) = 0, and satisfies the boundary condition D(r)(Pi) = 0 as required, but I haven't had any success with my first attempts in that direction.

 

Since the plot actually shows that the function a[1] is constant and equal to 0.1 you could plot the difference a[1](t)-0.1 as in this example where I use the logistic equation:
 

restart;
ode:=diff(x(t),t)=x(t)*(0.1-x(t));
res:=dsolve({ode,x(0)=0.1},numeric);
plots:-odeplot(res,[t,x(t)],0..10);
plots:-odeplot(res,[t,x(t)-0.1],0..10,caption="The difference between x(t) and 0.1");

The last plot:

I don't know exactly what was done in 3(d), but I would do the following, where the important part is to draw the orbit after it has had time to get close to the apparent attracting limit cycle. Thus I plot from t=200 to t=230:
 

restart;
with(plots):
EquationSet:= {diff(x(t),t) = -y(t) - z(t), diff(y(t),t) = x(t) + y(t)/5, diff(z(t),t) = 1/5 + (x(t) - 5/2)*z(t)};
Numsol1:=dsolve({EquationSet[],x(0)=0,y(0)=0,z(0)=0},[x(t),y(t),z(t)],numeric):
odeplot(Numsol1,[x(t),y(t),z(t)],200..230,numpoints=10000);
odeplot(Numsol1,[x(t),y(t)],200..230,numpoints=10000);
EquationSet2:= {diff(x(t),t) = -y(t) - z(t), diff(y(t),t) = x(t) + y(t)/5, diff(z(t),t) = 1/5 + (x(t) - 3)*z(t)};
Numsol2:=dsolve({EquationSet2[],x(0)=0,y(0)=0,z(0)=0},[x(t),y(t),z(t)],numeric):
odeplot(Numsol2,[x(t),y(t),z(t)],200..230,numpoints=10000);
odeplot(Numsol2,[x(t),y(t)],200..230,numpoints=10000);
## To produce an animation of the change from 2 to 3 via the value 5/2 I use the parameters option in dsolve/numeric.
##
EquationSetA:= {diff(x(t),t) = -y(t) - z(t), diff(y(t),t) = x(t) + y(t)/5, diff(z(t),t) = 1/5 + (x(t) - A)*z(t)};
res:=dsolve({EquationSetA[],x(0)=0,y(0)=0,z(0)=0},[x(t),y(t),z(t)],numeric,parameters=[A]);
Q:=proc(A) if not A::realcons then return 'procname(_passed)' end if;
      res(parameters=[A]);
      res
end proc;
animate(odeplot,[Q(A),[x(t),y(t),z(t)],200..230,numpoints=10000],A=2..3,frames=50,orientation=[0,20,30]);