@ramakrishnan Sure.

plots:-animate(DEplot,[ode,y(t),t=0..T,[seq([y(0)=y0,D(y)(0)=0],y0=-3..3)]],T=0.001..5);

Notice that I have T=0.001..5, not T=0..5. If you try the latter you will find that you get an error about the range being empty. This actually comes from DEplot. Compare
plots:-animate(plot,[sin(t),t=0..T],T=0..5); #OK
#and
DEplot(ode,y(t),t=0..0,[seq([y(0)=y0,D(y)(0)=0],y0=-3..3)]); #Errror

@BRUCELEE LinearAlgebra:-Dimension returns a sequence of two integers, not one. Try the following:

TransformationMatrix:=proc(CT::Matrix)
  local i,j,n,A;
  n:=LinearAlgebra:-Dimension(CT);
  print(n);
  seq(i,i=1..n);
 end proc:

TransformationMatrix(Matrix([[1,2],[3,4]]));

The solution is simple: Change Dimension to RowDimension.

@Kitonum Quite right. Initially I tried the inplace option as it has the advantage of not creating a new matrix A in each step. Here A defined as IdentityMatrix(N) won't work, but it does work when assigning to A as you are saying. I shall insert a note in my answer.
Thank you.

PS.  Your version without A is nice and fast (I tried setting n:=100).
This version with A defined without shape=identity is fast too:
for i from 1 to n do `.`(A,B[i],inplace) end do:

It avoids calling LinearAlgebra:-Multiply which probably is what slows this version down:
for i from 1 to n do LinearAlgebra:-Multiply(A,B[i],inplace) end do:

@ramakrishnan Yes, I see that in writing the response (in particular the line "Whether or not ...) I removed the assignment to res1 of the subs(_Z1=n,res) statement.
I have edited my answer above to include that assignment, but I may as well have done what you propose.

I cleaned up your lines and tried fsolve searching for nonnegative solutions, but it returned unevaluated.
I bring the code anyway, since it is easier to copy:

restart;
eqns := [A = (gr_c+delta)*kh^(1-alpha)/sav_rate, theta = Rk*(Rh-Rk)/(gamma*((Rh-Rk)^2+sigma^2)), theta = 1*1+kh, Rk = 1+rk-delta, Rh = 1+rh-delta, rk = A*alpha*kh^(alpha-1), rh = A*(1-alpha)*kh^alpha, sigma = sigmay/theta, varrho = Rap^((eps-1)*eps/(1-gamma)), Rap = Rk^(1-gamma)+(1-gamma)*Rk^(-gamma)*theta*(Rh-Rk)-.5*Rk^(-gamma-1)*gamma*(1-gamma)*theta^2*((Rh-Rk)^2+sigma^2), R = Rk+theta*(Rh-Rk), beta = ((1+gr_c)/R)^(1/eps)/varrho];
###
vals := [alpha = .36, delta = 0.6e-1, sigmay = sqrt(0.313e-1), gamma = 3, eps = .5, gr_c = 0.2e-1, sav_rate = .23];
###
eqns_a:=eval(eqns, vals);
vars:=indets(eqns_a,name);
fsolve(eqns_a, vars=~(0..infinity));

@sharena2 Do you have a reference for the statement that H is supposed to be monotonically decreasing?

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.