Here is an example:
 

restart;
ode:=diff(y(x),x)=y(x);
h:=0.1: #Stepsize
resRKF45:=dsolve({ode,y(0)=1},numeric,output=Array([seq(h*i,i=0..20)]));
resRK4:=dsolve({ode,y(0)=1},numeric,method=classical[rk4],stepsize=h, output=Array([seq(h*i,i=0..20)]));
A45:=resRKF45[2,1];
A4:=resRK4[2,1];
X:=A4[..,1];
Yd:=A4[..,2]-A45[..,2];
plot(X,Yd,caption=typeset("The difference between RK4 and RKF45\n with RK4 stepsize = ",h));

In this case you can also compare the two results with the exact solution y(x) = exp(x):
 

plots:-odeplot(resRKF45,[x,y(x)-exp(x)],0..2);
plots:-odeplot(resRK4,[x,y(x)-exp(x)],0..2);

@mrbayat There doesn't seem to be a difference between the mw and the mws file. Both have 2D-input, thus the mws-file isn't an mws-file.

But apart from that:
I executed your whole worksheet by using the !!! icon. Then I did:
 

indets(dsys,function);

What I saw includes weird constructions like X2(X2), diff(X2(X2), X2).

This must be cleared up.
The best is to remove X2 in N(X2) from the first line after restart so that you get:
 

Ws := (1/2)*N*k*T*(I1-3-2*log(J));

It appears as if a multiplication sign is missing in the line defining Wm. In any case it should be:

Wm := k*T*((J-1)*log(1-1/J)+chi-chi/J)/nu;

Finally, at the end you now get an error message saying that 7 boundary conditions are expected, but only two were provided. That is due to the fact that there are several (5) undefined constants in your system. Try

indets(dsys,name);

You get {L, T, X2, chi, k, nu}. So L, T, chi, k, nu needs to be given numerical values.

Have a look at PDEtools:-dchange.
For the present problem you can do as follows, where I have included initial conditions at t = t0.
First I'm doing ode and ics separately (more straightforward):
 

restart;
ode:=m*diff(x(t),t,t)+k*x(t)=0;
ics:=x(t0)=x0, D(x)(t0)=x1;
## Doing ode and ics separately:
PDEtools:-dchange({t=T*s, x(t)=X(s)*L},ode,[s,X]);
solve(%,{diff(X(s),s,s)});
subs(T^2=m/k,%);
ICS:=PDEtools:-dchange({t=T*s, x(t)=X(s)*L},{ics},[s,X]);
indets(ICS,function);
icsX:=subs(T=sqrt(m/k),solve(ICS,%));
##############
## Doing ode and ics at the same time:
res:=PDEtools:-dchange({t=T*s, x(t)=X(s)*L},{ode,ics},[s,X]);
## Separating the DE from the initial conditions
res1,res2:=selectremove(has,res,diff);
solve(res1,{diff(X(s),s,s)});
odeX:=op(subs(T^2=m/k,%));
indets(res2,function);
icsX:=subs(T=sqrt(m/k),solve(res2,%));

You may want to introduce s0 = t0*sqrt(k/m), X0 = x0/L, X1 = (x1/L)*sqrt(m/k).
While solving the new DE for the derivative(s) (here diff(X(s),s,s) ) is not necessary, it makes it easier to see (in general) which nondimensionalization to choose.
## You have introduced the constant L, which you can still choose to fit the situation. If x0 is never zero then L could be chosen as x0 so that X(s0) = X0 = 1.

Your second pade expression is when copied by me this:
 

pade*(theta, eta, [3, 3])

Obviously it should be:
 

pade(theta, eta, [3, 3])

With that change fsolve returns {b = -.26165340, d = .14500436}.

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.