@Carl Love Actually there is no complaint from dsolve/numeric when method=dverk78 and the ode is just

ode:= diff(y(x),x) = u; # No Re

The result is -0.555168769092963e-2+8.10096715234306*10^(-41)*I , which has the same real part as you get with Re(u) instead of u.
Methods ck45 and rosenbrock behave as rkf45.
No complaints from lsode or gear. No complaint from classical[rk4] (here without error tolerances of course).
### It appears that method=rkf45 doesn't complain when it can "see" from the start that the problem is complex (imaginary).
Examples:
 

restart;
ode:= diff(y(x),x) = 1+3*I; # A very visible I.
res:=dsolve({ode,y(0)=0},numeric,method=rkf45); #No problem
res(1); # Fine
##
g:=proc(x) if not x::numeric then 'procname(_passed)' else 1+3*I end if end proc;
odeg:=diff(y(x),x)=g(x); # I is not revealed at start
resg:=dsolve({odeg,y(0)=0},numeric,known=[g],method=rkf45); # Error
##
## In the following modified version I is visible from the start but is irrelevant if x<=10.
odeg2:=diff(y(x),x) = g(x) + piecewise(x>10,I,0); 
resg2:=dsolve({odeg2,y(0)=0},numeric,known=[g],method=rkf45); # OK
resg2(1); # OK
resg2(11); # OK

@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);