@vv Thank you for a very well chosen counter example to the "obvious" statement about the derivative of the trace in the paper:

Stability and bifurcation analysis in a predator–prey system with Michaelis–Menten type predator harvesting
by Dongpo Hu and Hongjun Cao.

@fatemeh1090 I understand that you had to remove the full pdf-file for copyright reasons.

Thus I won't be able to find the authors' reasoning.

I tried a little on my own, but didn't manage to get down to the result given for the derivative of the trace w.r.t. alpha.

I'll give the steps I tried:

restart;
local gamma;
##The system in the pdf-file. Right hand sides:
RHS:=[x*(1-x)-x*y,y*(delta-beta*y/x)-alpha*y/(gamma+y)];
ROsol:=solve(RHS=~0,{x,y});
RO:=op(indets(ROsol[2],specfunc(RootOf)));
E:=subs(ROsol[2],[x,y]);
da:=diff(E,alpha);
dxa,dya:=op(subs(RO=y,y=1-x,da));
J:=VectorCalculus:-Jacobian(RHS,[x,y]);
TR:=LinearAlgebra:-Trace(J);
TRxa:=normal(subs(y=1-x,x=x(alpha),TR));
diff(TRxa,alpha);
simplify(%,{TRxa=0}); # simplifying with the side relation TRxa=0
subs(diff(x(alpha),alpha)=dxa,x(alpha)=x,%);
res:=normal(%);

@fatemeh1090 The expression for the trace of J in (24) in the pdf-file is correct.

But how do you get from that to the expression given for the derivative of the trace w.r.t. alpha?
Of course x and y  depend on alpha, but it seems that the author is not using the explicit dependence we find easily with Maple.

The text continues after that with "Obviously, ..." and gives the derivative in terms of x etc.

Why is it obvious? There are several pages missing in the uploaded pdf-file. Maybe the derivative of x w.r.t. alpha has been found in those pages?

@vv I use a similar device for fsolve. I give the procedure (here PUV) the option remember to minimize the expense of calling the procedure (PUV)  several times.

@666 jvbasha Your worksheet BAHAM.mw doesn't run without several errors on my machine in Maple 2020, Maple 2018 and Maple 18 (which you apparently use).

Since this is a new project you should start a new thread, where everybody is invited to participate.

@666 jvbasha Please see my response above to your question directed to several users including me.

Besides that be aware of misunderstanding my answers about plotting errors.

In my answer comparing the solution of an ivp with the given bvp the initial conditions for the ivp (ics) were taken from the bvp. Thus the best I can get out of a comparison is the comfort that they are pretty close.
The ics are "contaminated" by the numerical errors in solving the bvp of course, so the ivp will be too, besides introducing its own errors. So we are not plotting the actual errors in solving the bvp.

@666 jvbasha Are you talking about iterations done by `dsolve/numeric/bvp` ?

If so the answer is that you have no control over the number of iterations done in order to satisfy the abserr requirement.

You can give a value to maxmesh in the range 32 ..13421772 according to the help. But that doesn't give you the control you (probably) want.

You need to write your own algorithm to gain that control.

@Kitonum Yes, it was introduced in Maple 2019.

@666 jvbasha Notice that in R_init T comes first, then diff(T(eta),eta), then f(eta), etc.

Thus use the order:
H, H1, F, F1, F2, G, G1, G2 := op(rhs~(R_init[2 .. -1]));

and everything else is fine.

@Alger What do you mean by saying that " Variables of the model should be null when C<40e-6 and can be not null when C>40e-6 with imposing high range." ?

Do you mean that for C <40e-6 all the dependent variables (like e.g. ids(t))  go to zero as t -> infinity?

And for C>40e-6 they tend to a periodic solution (not identically 0)?

@Alger I had another look at your system. Although it starts out looking linear it ends up being highly nonlinear. Thus there isn't any a priori reason to deny the actual existence of singularities.

The "singularity" problem I ran into with C = 0.000050 and range 0..50 disappeared, however, when I used your original equations eq1, eq2, eq3, and eq4 together with ep5 and ep6.

Unfortunately, for C = 0.000045 and range 0..100 I still got a "singularity" at the same point you got one.

There are obvious numerical difficulties with this problem, so I'm still not convinced that a genuine singularity is present.

Here is the revised worksheet:

MaplePrimes20-04-30_odesystemII.mw

@Alger When dsolve/numeric talks about a probable singularity one should (as you surely do) notice the adjective "probable".

What exactly is a singularity in this numerical context? What triggers the error (or warning) message?
Your system is very complicated, so it isn't easy to see. It would be nice to find a simpler example with the same problem.

@Alger You get a similar complaint when keeping C at 0.000050 but changing the range from 0..40 to 0..50:

Warning, cannot evaluate the solution further right of 47.288331, probably a singularity.

@arashghgood You ask in the title " Is there a solution?".

Clearly without additional requirements there are the very trivial solutions f(x,y,t) = F1(t)*y+F2(t) , (where F1 and F2 are arbitrary).

But I have no idea about how to solve this type of problem numerically.

You could set infolevel[int]:=5:

restart;
integrand1 := (x^3*exp(arcsin(x)))/sqrt(1 - x^2);
integrand2 := (x^3*exp(1)^arcsin(x))/sqrt(1 - x^2);
###
infolevel[int]:=5:
res1:=int(integrand1,x);

res2:=int(integrand2,x)

The first 4 steps are the same, but then they. Obviously the choice of approach is governed by the integrand and they are different.

exp(arcsin(x)) is of type function and exp(1)^arcsin(x) is of type `^` .