@Mariusz Iwaniuk I have no idea. If you are being sarcastic, peace be with you. I shall correct .2 to .1.
Thank you.

@MapleMathMatt You may be right. Mine is unchecked, so everything is allowed.

@dharr Since I do have a file at that place and readline works for me, let me give my output from
 

FileTools[Status]("F:/test.txt");

It is
[true, true, false, false, 1533641290, 68]

The two first show that I have read and write permission for that file according to the help for FileTools[Status].

@torabi Here is an animation that was made with the speedier version of fdsolve given in FracDiff.mw above:

The equations were taken from the paper in your link HOPFF.pdf.
The animation has 25 frames so it takes a while to make.
 

F:=omega*x-y^2;
G:=mu*(z-y);
H:=A*y-B*z+x*y;
params:={omega=-2.667,mu=10,A=27.3,B=1};
FGH:=unapply~(eval([F,G,H],params),t,x,y,z);
solfu:=alpha->if not alpha::numeric then 'procname(_passed)' else fdsolve(FGH,alpha,0..20,[.1,.1,.1],2000) end if;
plots:-animate(plots:-spacecurve,[solfu(alpha)[..,2..4] ],alpha=0.67..1);

@torabi Based on Rouben Rostamian's implementation I made a version that is considerably faster.
You can try it out. It appears to work fine, but I may work a little more on it at a later date.
I may have a look at your latest question later.

FracDiff.mw

@torabi q is just the value of N. Try

T := 20.0;
q := 1000;

and you get this plot:

I also plotted the 3D curve (corrected):

plots:-spacecurve([seq([sol[i][2],sol[i][3], sol[i][4]], i=0..q)],labels=[x,y,z]);

Actually I know nothing about fractional differential equations.

## Again it might be interesting or at least fun to compare with the corresponding ode system:
 

odesys:={diff(x(t),t)=F(t,x(t),y(t),z(t)),diff(y(t),t)=G(t,x(t),y(t),z(t)),diff(z(t),t)=H(t,x(t),y(t),z(t))};
ics:={x(0)=.39, y(0)=0.66,z(0)=-.4};
# Using params:
res:=dsolve(eval(odesys,params) union ics,numeric,maxfun=10^6);
plots:-odeplot(res,[x(t),y(t),z(t)],0..20,numpoints=1000);
plots:-odeplot(res,[x(t),y(t),z(t)],100..200,numpoints=10000);

@designay I have added the results in tabular form to my answer. Please have a look.

@Carl Love I have typesetting=standard and with that _EXPSEQ(1,-1) doesn't print as 1, -1. It just returns unevaluated.
 

restart;
interface(typesetting=standard);
interface(prettyprint) ; # 3
_EXPSEQ(1, -1); # returns unevaluated

Maple 2018.1, X86 64 WINDOWS, Jun 8 2018, Build ID 1321769
Standard Worksheet Interface, Maple 2018.1, Windows 10, June 8, 2018 Build ID 1321769

Sorry, but I cannot see anything in that image. Can you?
But don't bring images anyway. We need code either uploaded or presented here as text.

@umangbedi I tried letting pdsolve/numeric solve your system. It only allowed 4 of your 6 boundary conditions. Two of your pdes are first order in x and the third is second order.
 

restart;
###
A1 := 1.43727*10^9; A2 := 0.555556e-1; E := 80400; R := 8.314;
B1 := 2.61029; B2 := 0.555556e-1;
C1 := 159881; C2 := 2322.61; C3 := 1.0439*10^15;
###
PDE1 := A2*diff(C(x, t), t) = -diff(C(x, t), x)-A1*exp(-E/(R*Ts(x, t)))*(-0.2432e-1+1.07449*C(x, t)-0.21778e-1*C(x, t)^2-0.4266e-1*C(x, t)^3+0.1796e-1*C(x, t)^4);
PDE2 := B2*diff(Tg(x, t), t) = -diff(Tg(x, t), x)-B1*(Tg(x, t)-Ts(x, t));
PDE3 := C1*diff(Ts(x, t), t) = diff(Ts(x, t), x, x)+C3*exp(-E/(R*Ts(x, t)))*(-0.2432e-1+1.07449*C(x, t)-0.21778e-1*C(x, t)^2-0.4266e-1*C(x, t)^3+0.1796e-1*C(x, t)^4)+C2*(Tg(x, t)-Ts(x, t));
###                     
IC1 := C(x, 0) = 1;
IC2 := Tg(x, 0) = 900;
IC3 := Ts(x, 0) = 298; # You had Ts(0,t) (but used only the rhs, so OK).
bc2 := D[1](Ts)(0, t) = 0; bc1 := D[1](Ts)(1, t) = 0; # 2 for Ts
bc3 := Tg(0, t) = 900; bc4 := C(0, t) = 1; # 1 for Tg and 1 for C
### Not used in the pdsolve command here:
#bc5 := D[1](Tg)(1, t) = 0; bc6 := D[1](C)(1, t) = 0; 
res:=pdsolve({PDE1,PDE2,PDE3},{IC1,IC2,IC3,bc1,bc2,bc3,bc4},numeric,spacestep=0.02);
res:-plot3d(Ts,t=0..10);
res:-animate(Ts,t=10);

For your scheme it seems necessary to introduce the extra conditions bc5 and bc6 as is done here:
 

## Continuing without a restart:
xmin := 0;
xmax := 1;
tmax := 30;
N := 6;
##
dCdr := proc (h) (1/2)*(C[k+1](t)-C[k-1](t))/h end proc;
d2Tsdr2 := proc (h) (Ts[k+1](t)-2*Ts[k](t)+Ts[k-1](t))/h^2 end proc; 
dCdrb := proc (h) (C[k](t)-C[k-1](t))/h end proc; 
dTgdrb := proc (h) (Tg[k](t)-Tg[k-1](t))/h end proc; 
dTsdrb := proc (h) (Ts[k](t)-Ts[k-1](t))/h end proc; 
dTgdr := proc (h) (1/2)*(Tg[k+1](t)-Tg[k-1](t))/h end proc;
dTsdr := proc (h) (1/2)*(Ts[k+1](t)-Ts[k-1](t))/h end proc;
dCdrf := proc (h) (C[k+1](t)-C[k](t))/h end proc; 
dTgdrf := proc (h) (Tg[k+1](t)-Tg[k](t))/h end proc; 
dTsdrf := proc (h) (Ts[k+1](t)-Ts[k](t))/h end proc; 
##
h := xmax/(N-1);
##
stencil := subs(diff(Ts(x, t), x, x) = d2Tsdr2(h), diff(Tg(x, t), x) = dTgdr(h), diff(C(x, t), x) = dCdr(h), C(x, t) = C[k](t), Tg(x, t) = Tg[k](t), Ts(x, t) = Ts[k](t), x = k*h, {PDE1, PDE2, PDE3});
## Now introducing bc5 and bc4. 
bc5 := D[1](Tg)(1, t) = 0; bc6 := D[1](C)(1, t) = 0;
bcEqsExtra:=subs(k = N-1, {dTgdrb(h) = rhs(bc5),dCdrb(h) = rhs(bc6)});
bcEqs := subs(k = 0, {dCdrf(h) = rhs(bc4), dTgdrf(h) = rhs(bc3), dTsdrf(h) = rhs(bc2)}) union subs(k = N-1, {dTsdrb(h) = rhs(bc1)});
eqs := Vector(N-2): 
cnt := 0:
for k to N-2 do cnt := cnt+1; eqs(cnt) := stencil end do: 
##Removing curly brackets from each element in the vector eqs
eqs:=op~(eqs): 
eqs := [op(convert(eqs, list)), op(bcEqs),op(bcEqsExtra)]; 
ICs := []:
for k from 0 to N-1 do 
  ICs := [op(ICs), C[k](0) = rhs(IC1), Tg[k](0) = rhs(IC2), Ts[k](0) = rhs(IC3)] #Had Ts[0](t) instead of Ts[k](0)
end do; 
###
### An attempt:
resF:=dsolve([op(eqs), op(ICs)], numeric,abserr = 10^(-3), relerr = 10^(-3),stiff=true,output=Array([seq(tmax/100*i,i=0..100)]));
### Won't accept t values that high! 
### Trying t=0.84 instead of 30 (!)
resF:=dsolve([op(eqs), op(ICs)], numeric,abserr = 10^(-3), relerr = 10^(-3),stiff=true,output=Array([seq(0.84/100*i,i=0..100)]));
interface(rtablesize=20);
resF[1,1]; ## To see the meaning of the columns
A:=resF[2,1];
M:=Matrix(101,3*N,(i,j)->A[i,j+1]);
plots:-matrixplot(M[..,1..N]);
plots:-surfdata(M[..,1..N],0..0.84,0..xmax);
plots:-matrixplot(M[..,N+1..2*N]);
plots:-surfdata(M[..,N+1..2*N],0..0.84,0..xmax);
### All 3 at once:
plots:-display(Array([seq(plots:-matrixplot(M[..,k*N+1..(k+1)*N]),k=0..2)]));
plots:-display(Array([seq(plots:-surfdata(M[..,k*N+1..(k+1)*N],0..0.84,0..xmax),k=0..2)]));

So there are still some questions to be answered!

@tomleslie You wrote in your reply to Carl Love:
" I thought I was being somewhat more realistic than other "solutions" here because I didn't use the exact analytic solution to set up the shooting method. "

Well, in my answer, which predate yours, I only used the exact solution for comparison. I used a shooting method for the numerical solution. That I handled another ode (Problem 6) is irrelevant. Please have a look.

@tomleslie If you only have lines, points, and one filled area there doesn't seem to be covering in Maple 8 nor in Maple 2018.
But if you have two filled plots the one on top will cover the bottom one:
Define a smaller pink rectangle and display in two orders:
 

squPINK:=polygon([[1,1], [1,side], [side,side],[side,1]], color=pink):
plots[display]([squPINK,squ,poinC,tetA,lin[1],lin[2],lin[3],lin[4]], scaling=constrained,axes=normal);
plots[display]([squ,squPINK,poinC,tetA,lin[1],lin[2],lin[3],lin[4]], scaling=constrained,axes=normal);

In the first you see the pink rectangle, in the second you don't.
In Maple 8 you do see the black sides of the pink rectangle in the white plot though. Not in Maple 2018.

@9009134 I don't know if your system has a solution or not with the homogeneous BCS as given in my reply "Another look ...".
With BCS modified to BCS1 as below and using an approximate solution I made up a result turns up:
 

Digits:=15:
BCS1 := {theta(0) = 0, theta(1/1000) = 0, u(0) = 0, u(1/1000) = 0, w(0) = 0, w(1/1000) = 0, D(w)(0) = -0.1, D(w)(1/1000) = 0.1};
res:=dsolve(SYS union BCS1,numeric,approxsoln=[w(x)=100*x*(x-1/1000),u(x)=-3*10^(-7)*sin(2*Pi*1000*x),theta(x)=-200*(x-1/2000)],initmesh=2048,abserr=1e-7);

Plotting:
 

plots:-odeplot(res,[x,w(x)]);
plots:-odeplot(res,[x,u(x)]);
plots:-odeplot(res,[x,theta(x)]);

The first of these plots:

For the pure plot version I don't get red axes by the command:
 

plot(sin(x),x=-Pi..Pi,axis=[tickmarks=[color=red]], color=blue);
## Easier to see here:
plot(sin(x),x=-Pi..Pi,axis=[tickmarks=[color=red],thickness=3], color=blue);

I'm using Maple 2018.1 64bit with Windows 10.
Note: I can force the axes to be red (if I wanted to):

plot(sin(x),x=-Pi..Pi,axis=[tickmarks=[color=red],thickness=3,color=red], color=blue);

@9009134 I had another look at your system. I used the first ode to bring the total apparent differential order down from 10 to 8. That doesn't mean that I have a solution, but at least it clarifies the situation a bit:
 

restart;
dsys3 := {4.320502750*(diff(u(x), x, x))+4.320502750*(diff(w(x), x))*(diff(w(x), x, x)), 3.600418958*10^(-13)*(diff(theta(x), x, x))-1.206140351*theta(x)-1.206140351*(diff(w(x), x)), 1.206140351*(diff(theta(x), x))+1.206140351*(diff(w(x), x, x))+(4.320502750*(diff(u(x), x)+(1/2)*(diff(w(x), x))^2))*(diff(w(x), x, x)-8.379501384*10^20*(diff(w(x), x, x, x, x)))+(diff(w(x), x))*(4.320502750*(diff(u(x), x, x))+(diff(w(x), x))*(diff(w(x), x, x))-3.620365877*10^21*(diff(u(x), x, x, x, x))-1.086109763*10^22*(diff(w(x), x, x))*(diff(w(x), x, x, x))-3.620365877*10^21*(diff(w(x), x))*(diff(w(x), x, x, x, x)))+8.854000000*10^(-12)*(1.032500000-13000.00*w(x))/(2.500000000*10^(-6)-w(x))^2+6.503431892*10^(-32)/(2.500000000*10^(-6)-w(x))^4-4.451526315*10^10*(1.032500000-13000.00*w(x))*(diff(w(x), x))^2/(2.500000000*10^(-6)-w(x))^4+3.857989473*10^14*(diff(w(x), x))^2/(2.500000000*10^(-6)-w(x))^3-1.483842105*10^10*(1.032500000-13000.00*w(x))*(diff(w(x), x, x))/(2.500000000*10^(-6)-w(x))^3+9.644973683*10^13*(diff(w(x), x, x))/(2.500000000*10^(-6)-w(x))^2-1.089910330*10^(-9)*(diff(w(x), x))^2/(2.500000000*10^(-6)-w(x))^6-2.179820662*10^(-10)*(diff(w(x), x, x))/(2.500000000*10^(-6)-w(x))^5, theta(0) = 0, theta(0.1e-2) = 0, u(0) = 0, u(0.1e-2) = 0, w(0) = 0, w(0.1e-2) = 0, (D(u))(0) = 0, (D(u))(0.1e-2) = 0, (D(w))(0) = 0, (D(w))(0.1e-2) = 0}:
dsys3:=convert(dsys3,rational):
sys,bcs:=selectremove(has,dsys3,diff);
indets(sys,specfunc(diff));
sys[1];
eq:=solve(%,{diff(u(x),x,x)});
sys[3];
ode3:=eval(sys[3],eq);
indets(ode3,specfunc(diff));
sys[1..2];
SYS:={ode3,sys[1],sys[2]};
indets(SYS,specfunc(diff));
# Total order 2+2+4 = 8. So 8 boundary conditions, not 10.
bcs;
BCS:=remove(has,bcs,D(u)); # A natural choice 
dsolve(SYS union BCS,numeric); #No success
dsolve(SYS union BCS,numeric,method=bvp[midrich]); #No success