@mapleus Notice that although certainly a^n1*a^n2=a^(n1+n2) the identity you like to use is
a^(n1+n2) = a^n1*abs(a)^n2, which is not correct, but amounts to a reinterpretation of which integral you want to make zero. Also if n1 is even the problem is as before. Thus you should pick n1 odd to make sure the sign can change in the interval. But then a cleaner (but equivalent) interpretation for F(z) is
F(z)^n = signum(F(z))*abs(F(z))^n.

@Kitonum Thanks, you caught a special case: If h2=h1+1 you get a division by zero.
It seems that that can be handled by limits:

TresA:=limit(rhs(Tres2),h2=h1+1);
SresA:=limit(rhs(Sres2),h2=h1+1);
odetest({T(y)=TresA,sigma(y)=SresA},{Eq2,Eq3} union eval({bcs2},h2=h1+1));

@9009134 You can obtain the six solutions like this:

1. Pick as required {f2(0) = 0, f3(0) = 0, (D(f1))(0) = 0}. (Required by the differentiation resulting in ode).
Add to that one (only one) of {f2(1) = 0, f3(1) = 0, (D(f1))(1) = 0}. It doen't matter which since the remaining two will be satisfied because of conditions R1 and R2 (se above).
You now need to add 3 conditions from the available rest, i.e. those not having any of the conditions so far mentioned. Since there are 4 remaining conditions, you get 4 solutions that way.
2. Now still picking the required {f2(0) = 0, f3(0) = 0, (D(f1))(0) = 0} you add the 4 remaining conditions mentioned in point 1, i.e. this time you don't take any of {f2(1) = 0, f3(1) = 0, (D(f1))(1) = 0}. That gives you one solution. Thus so far you have 5 solutions.
3. Finally pick as required {f2(1) = 0, f3(1) = 0, (D(f1))(1) = 0}. If you pick any of {f2(0) = 0, f3(0) = 0, (D(f1))(0) = 0} you will have them all because of R1 and R2, so you won't be getting any new solutions. They will be the same as the 4 solutions from point 1.
So to get the last solution add to the required {f2(1) = 0, f3(1) = 0, (D(f1))(1) = 0} the 4 conditions left after removing {f2(0) = 0, f3(0) = 0, (D(f1))(0) = 0}.

You got 6 solutions.

@9009134 

Each of the following 13 sets

[{f1(0) = 0, f1(1) = 0, f2(1) = 0, (D(f3))(1) = 0}, {f1(1) = 0, f2(1) = 0, (D(f3))(0) = 0, (D(f3))(1) = 0}, {f1(0) = 0, f1(1) = 0, f2(1) = 0, (D(f3))(0) = 0}, {f1(0) = 0, (D(f1))(1) = 0, (D(f3))(0) = 0, (D(f3))(1) = 0}, {f1(1) = 0, (D(f1))(1) = 0, (D(f3))(0) = 0, (D(f3))(1) = 0}, {f1(0) = 0, f3(1) = 0, (D(f3))(0) = 0, (D(f3))(1) = 0}, {f1(0) = 0, f1(1) = 0, f3(1) = 0, (D(f3))(1) = 0}, {f1(0) = 0, f1(1) = 0, (D(f3))(0) = 0, (D(f3))(1) = 0}, {f1(0) = 0, f1(1) = 0, (D(f1))(1) = 0, (D(f3))(0) = 0}, {f1(0) = 0, f1(1) = 0, (D(f1))(1) = 0, (D(f3))(1) = 0}, {f1(0) = 0, f1(1) = 0, f3(1) = 0, (D(f3))(0) = 0}, {f1(1) = 0, f3(1) = 0, (D(f3))(0) = 0, (D(f3))(1) = 0}, {f1(0) = 0, f2(1) = 0, (D(f3))(0) = 0, (D(f3))(1) = 0}]

union

{f2(0) = 0, f3(0) = 0, (D(f1))(0) = 0}

and each of the following

[{f1(0) = 0, f2(0) = 0, (D(f3))(0) = 0, (D(f3))(1) = 0}, {f1(1) = 0, f2(0) = 0, (D(f3))(0) = 0, (D(f3))(1) = 0}, {f1(1) = 0, f3(0) = 0, (D(f3))(0) = 0, (D(f3))(1) = 0}, {f1(0) = 0, (D(f1))(0) = 0, (D(f3))(0) = 0, (D(f3))(1) = 0}, {f1(0) = 0, f3(0) = 0, (D(f3))(0) = 0, (D(f3))(1) = 0}, {f1(1) = 0, (D(f1))(0) = 0, (D(f3))(0) = 0, (D(f3))(1) = 0}, {f1(0) = 0, f1(1) = 0, f2(0) = 0, (D(f3))(0) = 0}, {f1(0) = 0, f1(1) = 0, (D(f1))(0) = 0, (D(f3))(1) = 0}, {f1(0) = 0, f1(1) = 0, (D(f3))(0) = 0, (D(f3))(1) = 0}, {f1(0) = 0, f1(1) = 0, (D(f1))(0) = 0, (D(f3))(0) = 0}, {f1(0) = 0, f1(1) = 0, f2(0) = 0, (D(f3))(1) = 0}, {f1(0) = 0, f1(1) = 0, f3(0) = 0, (D(f3))(1) = 0}, {f1(0) = 0, f1(1) = 0, f3(0) = 0, (D(f3))(0) = 0}]

union

{f2(1) = 0, f3(1) = 0, (D(f1))(1) = 0}

There actually only 6 different solutions. The reason is that these two relations:
R1:=-(D(f1))(0)+(D(f1))(1) = -(1/7)*f3(0)+(1/7)*f3(1)
R2:=-f2(0)+f2(1) = (11/7)*f3(0)-(11/7)*f3(1)
can found from
solve({sys[1],ode},{diff(f1(x),x,x),diff(f2(x),x)});
and then:
R1:=map(int,%[1],x=0..1,continuous);
R2:=map(int,%%[2],x=0..1,continuous);

Could you give us an explicit example?

To find solutions for which x>0,y>0 I believe it will be necessary to use numerical solution, but then of course you need all the parameters a, q1, ... to be concrete.

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