I wonder why you use `dsolve`instead of plain dsolve. There is no difference.

Doesn't the help page tell it as it is:

"The first integrals are set equal to generated global indexed variables K[i]  that denote arbitrary constants."

?

@tira I shall try again, but will begin with the conclusion:

Your expression V defined for y>0 has the limit 0 as y tends to zero from the right.

I do not use decimal numbers in order to show that the limit is indeed exactly zero.

restart:
t:=2/5; Pr:=71/100; Gamma:=1/5; Gr:=1/2;
## After this paste in your lines defining as you did all the terms v1, v2, ..., v8.
###
###Assuming that is done, continue with

C:=IntegrationTools:-Change(v4+v5-v3,u=y^2/x,[x]) assuming y>0; #Remember assuming
limit(C,y=0,right); #From the right in principle because the assumption y>0 was made
limit(v1-v2+C,y=0,right);
############Now we show that the 3 double integral terms v6, v7, v8 each has limit 0 as y->0:
## We shall use that
taylor(BesselI(1,x),x=0,4);
v6;
v6Int:=subs(int=Int,v6);
v6E:=simplify(eval(v6Int,BesselI=( ()->args[2]/2)));
# v6E is bounded above by
subs(y^2=0,v6E);
value(%);
v7;
v7Int:=subs(int=Int,-v7); #Notice minus: We are estimating absolute values
v7E:=simplify(eval(v7Int,BesselI=( ()->args[2]/2)));
# v7E is bounded above by
subs(y^2=0,v7E);
value(%);
v8;
v8Int:=subs(int=Int,v8);
v8E:=simplify(eval(v8Int,BesselI=( ()->args[2]/2)));
# v8E is bounded above by
subs(y^2=0,v8E);
value(%);
#Each of the upper bounds found for the absolute values of v6, v7, v8 has limit 0 as y->0.
##This finishes the proof of the conclusion stated in the beginning.

@tira The integrals in A-B of the same kind as v3 amount to v4+v5-v3.
Do the same substitution to that part of A-B :

IntegrationTools:-Change(v4+v5-v3,u=y^2/x,[x]) assuming y>0; ##Notice assumption
eval(%,y=0);
#Result 0.446573966
#  Then if you look at the contribution to A-B consisting of v1-v2, then you get
evalf(eval(v1-v2,y=0));
# except for roundoff likely the same.
#However, if y<0 is allowed in your setup, then you have a problem, because while the former result changes sign with y the latter does not.

Now left is the contribution to A-B from the double integrals v7+v8-v6. But because BesselI(1,x)  has a taylor expansion about 0 starting with x/2 each of the integrals itself is bounded in y thus with the factor y in front each of v6,v7,v8 tend to zero as y->0. 

Thus (except for the rounding error mentioned, which could be looked into) we see that the limit of A-B from the right as y->0 is 0.

The limit from the left is not zero though. But you may be assuming that y>=0 ?

@Kitonum It should be pointed out that your procedure Euler assumes that

1.  the derivative has been isolated on the left
2.  the dependent variable is called y
3.  the independent variable is called x

You have an h there, h=0.1. Does that mean that you wanted to solve the problem numerically?

I think the problem may be related to the fact that the right hand side of the equation for rho is not differentiable w.r.t. rho at rho=0.
Consider this simple example, where there is nonuniqueness when rho hits 0.

restart;
#Digits:=20: # Has no effect
ode:=diff(x(t),t)=-x(t)^(2/3);
res:=dsolve({ode,x(0)=1},numeric,events=[[x(t)-1e-14,halt]]);
plots:-odeplot(res,0..5); #Event triggered
res:=dsolve({ode,x(0)=1},numeric,events=[[x(t)-1e-15,halt]]);
plots:-odeplot(res,0..5); #Event not triggered, but singularity encountered.

@aureux Since I am showing both graphs in the same plot the wolf graph would be dwarfed by rhe elk graph. That was the only reason for multiplying by 100. You can obviously just remove that 100 or replace it with another number.
As an alternative you could make two graphs: One for elks and one for wolves.

@aureux
If your Maple version is not 18 then just remove the option size=[...] to odeplot. This option is new to Maple 18.
After that the plot should be fine.
You can stretch the plot with the mouse. This can be done by the size option in Maple 18.

Since the rest of your lines (starting with the assignment to eqs) refer to the nondimensionalized system and since you haven't defined that after the restart, they don't make any sense.

Notice that in my reply above I wrote:

"
################## Continuing with the dimensionless version ###########
However, if you would like to see how to get the same result from SYS in the previous exchange then continue from there with   "

where the bold face is in the original, i.e. I meant that it was important and should be noticed.

@THAPELO The replacement of N(t) with infinity in order to find the limits of X= <S,I1,I2,I3,I4> at infinity can be justified by noticing first that the resulting system for X is linear and autonomous.
It can be written

X'(t) = M.X+B  (*)

where M is a lower triangular matrix with negative numbers in the diagonal: Those are the eigenvalues of M.
B is a constant vector.

The actual system (keeping N(t)) can then be written as

X'(t) = M.X+B + F(t)  (**)

where F(t) -> zero vector as t->infinity.

This is enough for us to prove that the solutions of (**) tend to the equilibrium points of (*), which are simply found by
-M^(-1).B;
eval(%,params); #For the actual parameters.

The proof is basically the same as in the scalar case:

ode:=diff(x(t),t)=-a*x(t)+b+f(t), where a >0 and f(t)->0 as t->infinity

If need be I shall provide the proof for this scalar case, but it is rather straightforward.

@THAPELO As pointed out above there is no equilbrium with nonnegative coordinates for the system. Recall that an equilibrium solution is the same as a constant solution.
There is none even if we reformulate using AE.

Another question is: Does the solution [I1(t),I2(t),I3(t),I4(t),S(t),AE(t)] have a (finite) limit as t->infinity ?

The answer to that question I'm convinced is yes. The limit has I1,I2,I3,I4,S independent of initial conditions while the limiting value of AE depends on the initial conditions:

dsolve({Sys6E,AE(0)=AE0},{AE(t)});

The result has the form

Int(f(s)*exp(-a*s),s=0..t) +AE0;

where a=mu+nu and f = rho*(I1+I2+I3+I4).
The integral has a limit as t-> infinity, since f is bounded and a > 0.

@aureux Since I don't know your project and to what depth you are dipping into the questions of sustainability (a much overused term) I may have made the mistake of making this whole thing look more difficult than it is.

If you only want to see what actually happens for the parameters given and for the initial condition given then you need only do the following.

(I take for granted that when you wrote t=40 you meant t=0. The system is autonomous so the starting value for t may as well be 0. Anything else is just confusing.)

restart;
sys:=D(x)(t)=a*x(t)-a*x(t)^2/k-b*x(t)*y(t), D(y)(t) = -beta*y(t)+c*x(t)*y(t);
params:=[a = .71780, beta = .3067, b = 0.2072e-1, c = 0.9e-4, k = 17000];
sol:=dsolve({eval(sys,params),x(0)=6000,y(0)=40},numeric);
plots:-odeplot(sol,[[t,x(t)],[t,100*y(t)]],0..100,size=[1000,300]); #Notice 100*y
sol(100);

################## Continuing with the dimensionless version ###########
However, if you would like to see how to get the same result from SYS in the previous exchange then continue from there with

eqs:={X=k/(a*T),T=1/beta,Y=beta/b,k=(a/c)*A,a=B*beta};
params:=solve(eqs,{T,X,Y,A,B});

paramsConc:=[a = .71780, beta = .3067, b = 0.2072e-1, c = 0.9e-4, k = 17000];
ABTXY:=eval(params,paramsConc);
#Here I include initial values for u and v as parameters:
res:=dsolve({SYS,u(0)=u0,v(0)=v0},numeric,parameters=[u0,v0,A,B]);
#The parameters in the present concrete case:
L:=subs(ABTXY,[6000/X,40/Y,A,B]); #x(0)=6000, y(0)=40.
#Setting the parameters:
res(parameters=L);
#Readying input for odeplot so that it actually plots x and y although as functions of tau:
input:=subs(ABTXY,[ [[tau,X*u(tau)],[tau,100*Y*v(tau)]],tau=0..100/T]); #Notice 100*y
plots:-odeplot(res,op(input),size=[1000,300]);
###########Appearing as functions of t ######
##You can change the tickmarks and the labels to make the graphs look like graphs of x(t) and 100*y(t) like this:
input:=subs(ABTXY,[ [[tau,X*u(tau)],[tau,100*Y*v(tau)]],tau=0..100/T,
          tickmarks=[[seq(i*10/T=i*10,i=0..10)],default]]);
plots:-odeplot(res,op(input),labels=["t","x,100*y"],size=[1000,300]);

@THAPELO Noticing the exponential behavior of A and wanting to compute the values that I1,I2,I3,I4,S approach at infinity, I tried replacing A(t) by AE(t) = exp(-(mu+mu)*t)*A(t) for convenience.

Because of the long exchange, for the sake of clarity I paste the full code:

restart;
Digits:=15:
Sys1:=diff(S(t),t)=(1-h0)*Lambda1-(b*(I1(t)+phi1*I2(t)+phi2*I3(t))/N(t)+mu)*S(t);
Sys2:=diff(I1(t),t)=h1*Lambda1-(b*(I1(t)+phi1*I2(t)+phi2*I3(t))/N(t))*S(t)-(mu+rho+delta)*I1(t);
Sys3:=diff(I2(t),t)=h2* Lambda1+f*delta*I1(t)-(mu+theta+delta)*I2(t);
Sys4:=diff(I3(t),t)=h3*Lambda1+(1-f)*delta*I1(t)-(mu+theta+delta)*I3(t);
Sys5:=diff(I4(t),t)=h4*Lambda1+theta*(I2(t)+I3(t))-(mu+rho)*I4(t); ########* after theta
Sys6:=diff(A(t),t)=rho*(I1(t)+I2(t)+I3(t)+I4(t))+(mu+nu)*A(t); #minus? at the last term?
PDEtools:-dchange({A(t)=exp((mu+nu)*t)*AE(t)},Sys6,[AE]);
Sys6E:=combine(isolate(%,diff(AE(t),t)));
## Notice that the system ({Sys1,Sys2,Sys3,Sys4,Sys5,Sys6E} is no longer autonomous.
# Parameters
params:={Lambda1=0.029,b=0.5,phi1=0.25,phi2=1.01,mu=0.02,rho=0.1,delta=0.1,
  f=0.85,theta=0.2,nu=0.33,h0=0.4,h1=0.08,h2=0.02,h3=0.06,h4=0.4};
paramsE:=params union {subs(params,N(t)=S(t)+I1(t)+I2(t)+I3(t)+I4(t)+AE(t)*exp((mu+nu)*t))};
SYSE:=eval({Sys1,Sys2,Sys3,Sys4,Sys5,Sys6E},paramsE);
##Notice that AE(0)=A(0), so AE(0)=15.
Soln:=dsolve(SYSE union {S(0)=65,I1(0)=25,I2(0)=20,I3(0)=30,I4(0)=35,AE(0)=15},numeric);
vars:=[AE, I1, I2, I3, I4, S];
plots:-odeplot(Soln,[seq([t,v(t)],v=vars)],t=0..100,legend=vars,size=[1000,300],thickness=3);
#### We now try to find the values that AE, I1, I2, I3, I4, S seem to approach by other means than through Soln
#Since N(t)-> infinity as t -> infinity we do:
eval([Sys1,Sys2,Sys3,Sys4,Sys5],N(t)=infinity);
evalindets(%,function,s->op(0,s));
solve(%,vars[2..6]);
eval(%,params);
subs(Soln(500),vars[2..6]=~vars[2..6](t));
## Notice the perfect agreement except for S which is close though.
## The reason is that S(t) decreases rather slowly towards its value at infinity, as exp(-mu*):
#We look at that
eval(Sys1,N(t)=infinity);
dsolve(%);
eval(%,params);
plots:-odeplot(Soln,[t,ln(S(t)-0.87)],t=100..200);
#The slope ought to be -mu = -1/50 which appears correct.

@THAPELO Please notice that I edited the reply above.

@THAPELO Substituting parameters first helps make the two solutions look much less complicated:

Although not necessary for solve I use the opportunity to replace A(t) by A, I1(t) by I1, etc. by the use of evalindets:

eqs:=evalindets( eval(rhs~({Sys1,Sys2,Sys3,Sys4,Sys5,Sys6}),params),function,x->op(0,x));
solve(eqs,{A,S,I1,I2,I3,I4});

pts:=map(subs,[%],[A,S,I1,I2,I3,I4]); #lower case m in map.

You got 2 equilibrium points none of which consists entirely of nonnegative values. Irrelevant?

However, notice that A(t) grows exponentially with t as exp((mu+nu)*t) in fact.

But the other functions appear to approach positive values as t -> infinity. Try plotting those only (i.e. don't include A(t)).