@9009134 I'm assuming that your system is a mathematical model of some physical setup. Thus this would determine the boundary conditions.
Notice that the system is first order in f2, thus you cannot have these: D(f2)(1) = 0,  D(f2)(0)=0.
This leaves 10 on your list.
From a purely mathematical point any set of 7 boundary conditions chosen from the 10 you write down are worth a try. You need to keep in mind though that as stated earlier either
10*f2(0)+12*D(f1)(0)+14*f3(0) = 0
or
10*f2(1)+12*D(f1)(1)+14*f3(1) = 0
has to be satisfied

As an example you could pick 2 conditions for f1, 1 condition for f2, and 3 conditions for f3.
Concrete example, which will satisfy 10*f2(0)+12*D(f1)(0)+14*f3(0) = 0:
bcs1:=f1(1) = 0,D(f1)(0) = 0 ,f2(0) = 0 ,f3(1) =0,f3(0) =0 , D(f3)(1) = 0,  D(f3)(0)=0;

Notice though that it is not necessary that the number of conditions for a particular unknown (f1,f2, or f3) match its highest order. A very simple example illustrates that:
restart;
sys:=diff(y(x),x,x)-y(x)=0, diff(z(x),x)=y(x);
dsolve({sys}); #General solution
dsolve({sys,y(0)=0,z(0)=1,z(1)=0}); #Exact solution of bvp
dsolve({sys,y(0)=0,z(0)=1,z(1)=0},numeric); #Numeric solution of bvp

@ramakrishnan

1. DEplot is a procedure in the DEtools package. If you want to use its short form: DEplot you need to load the package:
with(DEtools):
if you prefer you can use the long form instead:  DEtools[DEplot].
2. The animate command must use dl and not ode, since dl is the name you gave the equation.

Yes it is possible to upload a worksheet. Use the thick green arrow icon in the editor.
Another way to present reasonably readable code is to remove all output from the worksheet before copying. I almost always do that myself: In Maple, go to Edit/Remove Output/From Worksheet.

@ramakrishnan I suggest you show the attempts you made even if they failed. Then we might be able to get what you want.

@Markiyan Hirnyk I think we have had this discussion before. But I repeat what I must have stated at that time: the constant coefficient case is straightforward.

@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