@mapleus The problem is that M_F is negative between 0 and 2. This causes no problem when n is an integer, but if is a float then Maple uses the usual definition of a^n:

a^n = exp(n*ln(a) )

The (complex-) valued logarithm in Maple is the principal branch, i.e.

ln(a) = ln(abs(a)) + I*arg(a), where  -Pi < arg(a) <=Pi
and the ln on the right is the usual real-valued logarithm.

Thus (-1)^4.7 = exp(4.7*ln(-1)) = exp(4.7*I*Pi) = -0.587785+0.809017*I.

To compensate you need to start with large (positive) values of X_1 and X_2 since M_1 and M_2 are positive.
Certainly X_1 = 1000 and X_2 = 1000 makes fsolve start and go for a while. Whether you ever get success I don't know.

However, since M_1 and M_2 are nonnegative on the interval 0..2 the quantity (M_F+X_1*M_1+X_2*M_2)^n cannot be nonnegative on 0..2 if one_int and two_int are to be zero.
Thus values for X_1 and X_2 that makes M_F+X_1*M_1+X_2*M_2 > 0 for all z in 0..2 will not solve your problem.
If you don't reformulate your problem you cannot avoid the appearance of imaginary numbers.

You have two integrals of the form
A=Int(F(z)^n*h(z),z=0..2)
where F(z) is real and h(z)>=0 for all z. n is (say) 4.7.
The imaginary part of the integral is
Ai = Int( abs(F(z))^n*sin(n*arg(F(z)))*h(z),z=0..2);
On that part of the interval in which F(z) > 0 we have arg(F(z)) = 0. Hence the integral over that part is zero.
On the part where F(z) < 0 we have arg(F(z)) = Pi. Hence the integral over that part is nonzero since sin(n*Pi)<>0 unless F(z)*h(z) = 0 on that part.
The real part of the integral is
Ar = Int( abs(F(z))^n*cos(n*arg(F(z)))*h(z),z=0..2);
On the integral over that part where F(z) > 0 is nonzero since arg(F(z))=0 and cos(0)=1 unless F(z)*h(z)=0 on that part.
Combining the two results: The integral A can only be zero if F(z)*h(z) = 0 for (almost) all z.

@Kitonum Starting with an expression with indexed variables:
n:=5:
u:=add(A[i]*x^i,i=0..n);
u1:=evalindets(u,indexed,s->cat(op(0,s),op(s)));

@mapleus In Maple it is often so that if an equation is expected (or accepted), but an algebraic expression u is given as input, then this is taken to mean u=0.
Simple example:
solve(x+1,x);

With procedural input to fsolve it is similar except here equations are not accepted.
Thus if you wanted evalf(Int(B*(M_F+X_1*M_1+X_2*M_2)^n*M_1,z=0..2*l)) = 5 you would change the procedure to

one_int:=proc(X_1,X_2) local res;
  res:=evalf(Int(B*(M_F+X_1*M_1+X_2*M_2)^n*M_1,z=0..2*l));
  print([_passed],res);
  res-5  #5 subtracted
end proc;

The same should be done to two_int.

@mapleus 
1. For each pair of concrete values for X_1 and X_2 one_int and two_int each returns a number. Thus one_int and two_int are functions of X_1 and X_2. z is just an integration variable and cannot be an input to one_ int and two_int.
2. The combination evalf(Int(....)) evaluates the integral numerically. Lower case 'int' attempts exact integration, which I believe has no chance at all in this case.
3. Surely the help for fsolve gives the syntax for procedural input.
Look in particular under ?fsolve,details
The first argument here is a list of the two procedures one_int and two_int.
The second argument I used (optional) is a list of two starting values (guesses). Here intervals could have been given instead as in [-2..0,0..1] or whatever is relevant.

@9009134 If you insist on satisfying 10*f2(0)+12*D(f1)(0)+14*f3(0) = 0 and by using only conditions from
bcs:=f1(1) = 0,f1(0) = 0 , D(f1)(1) = 0,  D(f1)(0)=0,f2(1) = 0,f2(0) = 0 ,f3(1) =0,f3(0) =0 , D(f3)(1) = 0,  D(f3)(0)=0;
then we have the requirement
bcs0r:={D(f1)(0)=0,f2(0) = 0 ,f3(0) =0};
Now you have to pick 4 conditions from the remaining 7:
bcs0a:={bcs} minus bcs0r;
There are binomial(7,4)=35 different choices.
I went though these and found immediate results for 13 of these choices. By comparing those 13 it turned out that there were only 5 different ones.
##Edited:
By doing the similar thing for 10*f2(1)+12*D(f1)(1)+14*f3(1) = 0 out of the 35 choices again 13 gave immediate results. But only one new solution was found. Thus altogether 6 different solutions were found.
When I say "immediate result" I mean that I automated this and let the process go on to the next bvp if dsolve returned an error.

@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?