@adel-00 eq1 is a polynomial equation of degree 3. Maple gives 3 solutions when it realizes that fact and when they are real. That can complicate matters.
Here I find all three (complex) solutions for u and after that for v.

restart;
eq1:=(alpha+(l+alpha)*u+alpha*k*u^2)*a=
u*(alpha+(l+alpha)*u+alpha*k*u^2)*(1+l*alpha*b/((alpha+(l+alpha)*u+alpha*k*u^2)));
eq2:=a=(u*(1+l*v/(1+u+k*u^2)));
factor((rhs-lhs)(eq1));
eq1:=collect(%,u); #Now obiously a polynomial
params:={l=1,alpha=1,b=100,k=50};
#Finding 3 solutions for each a. I have made the increment in a 0.1 to see the changes more clearly.
cupoints:=[seq(fsolve(eval(eq1,params),u,complex),a=0..20,.1)]:
plots:-complexplot(cupoints,style=point);
#Making an animation:
Su:=seq(plots:-complexplot(cupoints[1..n],style=point,symbolsize=15),n=1..nops(cupoints)):
plots:-display(Su,insequence=true);
#Turning to the v-values (slightly repetitive):
cpoints:=[seq([[a,a,a],[fsolve(eval(eq1,params),u,complex)]],a=0..20,.1)]:
Vp:=eval(solve(eq2,v),params);
Sau:=[seq( zip~( (x,y)->[a=x,u=y],op(cpoints[i])),i=1..nops(cpoints))]:
cvpointsTemp:=ListTools:-Flatten([seq(eval~(Vp,Sau[i]),i=1..nops(cpoints))]):
#Removing disturbing Float(undefined):
cvpoints:=remove(type,cvpointsTemp,undefined):
plots:-complexplot(cvpoints,style=point);
#Animating the v-points
Sv:=seq(plots:-complexplot(cvpoints[1..n],style=point,symbolsize=15),n=1..nops(cvpoints)):
plots:-display(Sv,insequence=true);

Added: A simple graphical illustration showing when there are 1, 2, or 3 real solutions:

plots:-animate(plot,[(eval(eq1,params)),u=-1..7,-20..20],a=2.7..6.5);

@adel-00 a=0 is a little special, so I left that out:

vpoints:=seq([a,subs(fsolve(eval({eq1,eq2},params),{u,v}),v)],a=1..20);
plot([vpoints]);
plot([vpoints],style=point);

# a = 0 is special because then the equations are satisfied if u = 0 and v is any number:

eval({eq1,eq2},params union {a=0,u=0});

#A slightly different approach:

V:=aa->subs(fsolve(eval({eq1,eq2},params union {a=aa}),{u,v}),v);
plot(V,0..20);

@adel-00 a=0 is a little special, so I left that out:

vpoints:=seq([a,subs(fsolve(eval({eq1,eq2},params),{u,v}),v)],a=1..20);
plot([vpoints]);
plot([vpoints],style=point);

# a = 0 is special because then the equations are satisfied if u = 0 and v is any number:

eval({eq1,eq2},params union {a=0,u=0});

#A slightly different approach:

V:=aa->subs(fsolve(eval({eq1,eq2},params union {a=aa}),{u,v}),v);
plot(V,0..20);

I tried this before I saw yours, but my version fails:

f:=proc(x,y) if not type([x,y],list(numeric)) then return 'procname(_passed)' end if;
      `if`(y>0,[x,4*y],[x,y])
end proc;

Edited: In the initial version of this comment I forgot 'return' . But the edited is the one intended and it doesn't work.

Edited once again: The version works if  _passed is replaced by x,y. Probably because _passed will be arguments passed to another procedure.

I tried this before I saw yours, but my version fails:

f:=proc(x,y) if not type([x,y],list(numeric)) then return 'procname(_passed)' end if;
      `if`(y>0,[x,4*y],[x,y])
end proc;

Edited: In the initial version of this comment I forgot 'return' . But the edited is the one intended and it doesn't work.

Edited once again: The version works if  _passed is replaced by x,y. Probably because _passed will be arguments passed to another procedure.

You say that you would like to find the analytic expression for that sum.

It may be that none such exists!

I haven't looked at this version, but have just answered your additional question in your previous posting. Please take a look at that.

@felixiao I made some changes. There were still some instances of F(tau), P(tau), Q(tau)  which needed to be replaced with F, P, Q, respectively. I replaced 'matrix' with 'Matrix' and `&*` with `.`. Also 'Inverse' was replaced with `^(-1)`.

MaplePrimes13-03-08.mw

@felixiao I made some changes. There were still some instances of F(tau), P(tau), Q(tau)  which needed to be replaced with F, P, Q, respectively. I replaced 'matrix' with 'Matrix' and `&*` with `.`. Also 'Inverse' was replaced with `^(-1)`.

MaplePrimes13-03-08.mw

I was briefly looking at this and decided to remove all the underscores. That should be harmless of course.
Then tested my answer to your previous question, but now I got 94 solutions instead of 90. Weird!
I used Edit/ Find/Replace in Maple.  For the record, here is the code with the underscores removed:
restart;
S := Matrix( [ [a1, a2, a3, a4], [b1, b2, b3, b4], [c1, c2, c3, c4], [d1, d2, d3, d4] ] );
eq1:=a1*d1 = b1*c1, a2*d2 = b2*c2, a3*d3 = b3*c3, a4*d4 = b4*c4;
R := Matrix( [ [s1, t1, r1, l1], [s2, t2, r2, l2], [s3, t3, r3, l3], [s4, t4, r4, l4] ] );
eq2:=s1*l1 = t1*r1,s2*l2 = t2*r2,s3*l3 = t3*r3,s4*l4 = t4*r4;
T := Matrix( [ [1, 0, 0, 1], [0, 1, 0, 0], [0, 1, 0, 0], [1, 0, 0, 1] ] );
S.R=~T;
convert(%,set);
res:=solve(% union {eq1,eq2});
nops([res]);
res[1];

Could you give us the expression you want to minimize? In text here or as an uploaded worksheet.

restart;
sys := [diff(y(x), x) = -(4*cos(x)*y(x)+z(x)*cos(x)^2+3*z(x))/(sin(x)*(cos(x)^2-9)),
diff(z(x), x) = -(y(x)*cos(x)^2+3*y(x)+4*z(x)*cos(x))/(sin(x)*(cos(x)^2-9))];

normal(sys[1]+sys[2]);
#Letting u = y+z:
eqS:=eval(%,{z=0,y=u});
U:=dsolve(eqS) assuming real;
normal(sys[1]-sys[2]);
#Letting v = y-z:
eqD:=eval(%,{z=0,y=v});
dsolve(eqD) assuming real;
V:=subs(_C1=_C2,%);
y(x)=rhs((U+V)/2);
z(x)=rhs((U-V)/2);

@Markiyan Hirnyk I noticed that the right hand sides of sys shared singularities defined by sin(x) = 0.
Thus the variable change from x to t given by diff(x(t),t)=sin(x(t)) seemed natural because the ode system in X(t), Y(t) would be without singularities (by the chain rule dy/dt = (dy/dx)*(dx/dt)=(dy/dx)*sin(x) ).
Omitting the last line in my comment above and setting
res:= collect(%,[_C1,_C2],normal);
#the solutions can be extended to all of R by

res2:=eval(res,{(1+cos(x))^(-1/2)=1/(sqrt(2)*cos(x/2)),sqrt(1+cos(x))=sqrt(2)*cos(x/2)});
odetest(res2,sys); #Checking
y1,z1:=op(eval(subs(res2,[y(x),z(x)]),{_C1=1,_C2=0}));
y2,z2:=op(eval(subs(res2,[y(x),z(x)]),{_C1=0,_C2=1}));
plot([y1,y2,z1,z2],x=0..4*Pi);

Since you left sol1[2] to the MaplePrimes users I thought I might try a variant of the idea.

restart;
sys := [diff(y(x), x) = -(4*cos(x)*y(x)+z(x)*cos(x)^2+3*z(x))/(sin(x)*(cos(x)^2-9)),
diff(z(x), x) = -(y(x)*cos(x)^2+3*y(x)+4*z(x)*cos(x))/(sin(x)*(cos(x)^2-9))];
#Making the change of variable given by sin(x)*diff(y,x) = diff(y,t),
# i.e. diff(x(t),t)=sin(x(t)) and choosing x(0) = Pi/2:  
int(1/sin(x1),x1=Pi/2..x, AllSolutions) assuming x>0,x<Pi;
St:=t=combine(%) assuming x>0,x<Pi;
dsolve({diff(x(t),t)=sin(x(t)),x(0)=Pi/2});
Sx:=subs(x(t)=x,%);
PDEtools:-dchange({Sx,y(x)=Y(t),z(x)=Z(t)},sys,[t,Y,Z]):
sys2:=combine(solve(%,{diff(Y(t),t),diff(Z(t),t)}));
tres:=dsolve(sys2) assuming t>0;
PDEtools:-dchange({St,Y(t)=y(x),Z(t)=z(x)},tres,[x,y,z]):
simplify(%) assuming real;
collect(%,[_C1,_C2],normal);
subs(sin(x)=sqrt(1-cos(x))*sqrt(1+cos(x)),%);

@matja There must have been something else interfering with the assume command since the following works:

restart;
assume(a>0);
limit(exp(-a*x), x = infinity);