I had no problem with Gr=7, but had to work to get results for Gr=7.1,..,7.7. No immediate success after that.
Insert: with method=bvp[midrich] I could go to Gr=7.6 without any problem or hard work.
Here is an animation showing F', G, H for those Gr-values (colors: red,blue,green).




Notice that H (green) changes a lot. In particular D(H)(0) increases rapidly. If there is a real mathematical problem here or just a numerical one I have no idea.

@9009134 The system sys2 is of total order 2+1+4 = 7, thus you can and must have 7 boundary conditions. These boundary conditions can involve f1, D(f1), f2, f3, D(f3), (D@@2)(f3), (D@@3)(f3) at the two boundary points x=0 and x=1.
In order for sys2 to have the same solution (f1,f2,f3) as sys the boundary conditions must include condition(s) sufficient to make either
10*f2(0)+12*D(f1)(0)+14*f3(0) = 0
or
10*f2(1)+12*D(f1)(1)+14*f3(1) = 0
satisfied.
One of these could be included as it is stated here, and then you need 6 other independent conditions.
If 3 of your conditions are f2(1)=0, D(f1)(1) =0, f3(1) = 0 then you should not add any of the above two conditions.

@9009134 In view of my previous comment and further exploration I think you are demanding too much in bcs. In other words, I believe that your original problem has no solution.
Suppose you keep the requirements f2(1) = 0 and f3(1) = 0. Then it follows from sys[2] that D(f1)(1)=0. Thus that must be one of the conditions. That makes it necessary to remove some other condition. The condition removed will not be satisfied.
Examples:
#Leaving out D(f3)(1)=0, but adding D(f1)(1)=0:
bcs2:={f1(0) = 0, f1(1) = 0, D(f1)(1)=0, f2(1) = 0, f3(0) = 0, f3(1) = 0, D(f3)(0) = 0};
res3:=dsolve(sys2 union bcs2,numeric,output=listprocedure);
plots:-odeplot(res3,[[x,100*f1(x)],[x,f2(x)],[x,f3(x)]],0..1);
res3(1);
#Leaving out D(f3)(0)=0, but adding D(f1)(1)=0:
bcs2a:={f1(0) = 0, f1(1) = 0, D(f1)(1)=0, f2(1) = 0, f3(0) = 0, f3(1) = 0, D(f3)(1) = 0};
res4:=dsolve(sys2 union bcs2a,numeric,output=listprocedure);
plots:-odeplot(res4,[[x,100*f1(x)],[x,f2(x)],[x,f3(x)]],0..1);
res4(0);
#Leaving out f3(0)=0, but adding D(f1)(1)=0:
bcs2b:={f1(0) = 0, f1(1) = 0, D(f1)(1)=0, f2(1) = 0, f3(1) = 0, D(f3)(0)=0,D(f3)(1) = 0};
res5:=dsolve(sys2 union bcs2b,numeric,output=listprocedure);
plots:-odeplot(res5,[[x,100*f1(x)],[x,f2(x)],[x,f3(x)]],0..1);
res5(0);

@9009134 You have a good point here.
Supposing that the original system sys with bcs has a solution for f1,f2,f3. Then certainly that solution satisfies the new system sys2 still with bcs.
The problem is going the other way: We have a solution to sys2 with bcs. Does this solution satisfy sys, which here means: does it satisfy sys[2]?
sys[2] is 10*f2(x)+12*(diff(f1(x), x))+14*f3(x) = 0. From ode it follows by integration that
10*f2(x)+12*(diff(f1(x), x))+14*f3(x) = C for some constant C. From bcs we don't have enough information to conclude that the value of C = 0, indeed it might not be.
An easy check reveals that in fact C = 0.005021010301197:
plots:-odeplot(res,[x,lhs(sys[2])],0..1);

Thus we must conclude that if sys with bcs has a solution, then sys2 with bcs has at least 2 solutions.

@9009134 This is what I did:

MaplePrimes14-11-08odesysINT.mw

@9009134 I have no problem executing your statements till and including the ones involving PolynomialInterPolation.
I left out the solve statements since their results are not assigned to, and not used. When I had similar solve statements in the link given above, they where there for the purpose of seeing if solution for the highest derivatives was possible. That was not the case before differentiation, but was the case after.
However, what follows after the PolynomialInterpolation must be missing something or is left over from your previous attempts. There I see a reference to dsys3, which is never defined in that worksheet.

The plots
plots:-odeplot(res,[[x,f1(x)],[x,f2(x)],[x,f3(x)]],0..1);
and
plot([fy11,fy22,fy33],x=0..1);

are (not surprisingly) very similar (indistinguishable). I suppose that your reason for doing the interpolation is that you want analytical expressions for the approximate solutions.

Just from picking rather arbitrary values for the parameters you get
eval(cdm_ode,{c0=.2, n=5,sigma=175, s0=200, ks=.5, h1=1000, h2=.1,kp=0.5, A=1, B=2});
#It is not very surprising then that
err(.2,5,200,.5,1000,.1,.5,1,2);
detects a singularity.
Trying the presumably simplest case n=1:

sys:=eval(cdm_ode,{c0=.2, n=1,sigma=175, s0=200, ks=.5, h1=1000, h2=.1,kp=0.5, A=1, B=2});
res := dsolve([sys, y1(0) = 0, y2(0) = 0, y3(0) = 0, y4(0) = 0, y5(0) = 0, y6(0) = 175], numeric, range = 0 .. tol_t);
plots:-odeplot(res,[t,y1(t)],0..10);
#The singularity seems obvious.

@Dmitry You cannot know that it always has an additional solution, but you cannot know that it doesn't either.
eq1:=simplify(eval(eq,params=~{seq(1..nops(params))}));
res:=solve(eq1,t3);
evalf[50](eval(res,t1=1));
plots:-implicitplot(eq1,t1=0..5,t3=0..5,gridrefine=3,grid=[50,50]);
plot(eval(lhs(eq1),t1=1),t3=0..2,-1..1);
Digits:=20:
fsolve(eval(eq1,t1=1),t3=1.33);

#So here there are clearly 2 positive solutions.

@Dmitry No, on the contrary. I mean that generally you should expect more than one solution (even more than one real solution). I added a few lines above which shows that when all parameters are 7 and t1=1 then besides t3=1 and t3 = 1.067641139 there is also a complex solution t3 = -1.238855445-.8386023959*I, and even one more real solution.

You could try following the advice given, i.e. using optional arguments
maxmesh=256 (or higher, 8192 being the highest available) and/or abserr=1e-3 (just an example).
That may or may not help. If not, let us see the whole problem.

@Vesnog Here is a simple example:

eq:=tan(x)=a*x;
sol:=[seq(fsolve(eq,x=0..Pi/2,avoid={x=0}),a=[1.1,1.6,Pi])];

Your expression RHS-LHS is not a polynomial in x, thus fsolve will not return all 3 solutions on Pi/2..lim.
The method of splitting you are using works fine, but obiously has the disadvantage of relying on the graph or at least on knowledge of the singularities of RHS-LHS.

@edrobina I have no problem with your worksheet in Maple 18.

So apparently the set  {a, b, c, d, e, f, g, h, i, j, k, l, m, n, o, p} does not only have (unassigned) names a,b,c,...
What is actually in {a, b,  c, d, e, f, g, h, i, j, k, l, m, n, o, p} ?

I think this was a display prpblem in Maple 16 (later corrected as Markiyan points out).
To attempt a verification of this you may try:
taylor(1/(1+z^2), z = 1, 4);
convert(%,polynom);

That it is a display problem is indicated by the fact that the term 1/2 and the term -1 in the incorrect looking result would automatically be simplified to -1/2.