@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

Instead of an image could you upload a worksheet? Use the big green arrow in the MaplePrimes editor.

@samen You may want to compare your results to what you get by applying Maple's pdsolve/numeric.
So what is the pde?
Example which may be close to yours.
 

restart;
pde:=diff(v(x,t),t)=a*diff(v(x,t),x,x);
res:=pdsolve(eval(pde,a=10^(-6)),{v(x,0)=1+x*0.9,v(0,t)=0.1,v(1,t)=0.5},numeric,timestep=.0001,spacestep=0.1);
res:-plot3d(t=1);

You can also give a method (including CrankNicholson). See the help page
?pdsolve,numeric,method

Note: In your code you should change the definition of the matrix nn to
 

nn := Matrix(N+2, N+2,(i, j)-> u[i-1, j-1]):

otherwise you won't see the boundary at x=1.
 

@9009134 I think that the immediate problem of not being able to convert to an explicit first order system can be overcome if you are willing to change the boundary conditions somewhat.

But you also have a possibly serious problem in the interior of your interval 0 .. 1/1000.
The denominators (1/400000 - w(x) )^n, n= 2..6 may easily become zero since w(0)=0 and w(1/1000) =0.