@Carl Love 

Using fsolve(..., complex) as you suggested gives a solution close to yours to within  10^-11 (or less) which verifies rel to within 10^-42.
Given that diff(Re(rel[1]), A) is infinite at this point, the slightest variation in A makes "huge" differences in the values of rel[1].

2023_vs_2015.mw

So you're right, there is indeed a second real solution

@Carl Love 

Thanks Carl... but this is not the result I get with Maple 2015:

eval(rel, {
    A = -2.7553365135418814642586082436429575890825402826031,
    B = -0.70285804987973303586180028708027467941012949957141
})

[                                                     -43    
[0.00004788232651393033381187767465396229938775 - 2 10    I, 

                                                        -40
  1.8649856419668477410903757547200476549949700479584 10   

                                                           -41  ]
   - 1.2309622539151182340327668587211822233603338843467 10    I]

A version issue?

@NIMA112 


Are there other real roots than the one @Rouben Rostamian got?
I don't think so (the couple (A, B) @Carl Love found doesn't verify the equations with an enough small error to be considered, IMO, as a solution).

Some details are given here Fung_sand15.mw

If the target starts from the left focus with velocity V_T and the predator from the left vertex with velocity V_P then: the predator will catch the prey only if V_P >V_T ant it will catch it at the center of the ellipse iif V_P / V_T = 1/e (e=eccentricity, assumed > 0 and < 1).
This is a particular situation where the capture at the center of the ellipse is not the rule.

So I doubt that, without extra conditions, the capture always takes place at the center of the ellipse

@RezaZanjirani 

Suppose you have to functions f[1](x) and f[2](x).
Heaviside(f(x)) (let's put aside the "x=0 case") has two values: 0 if f[1](x) < 0 ans 1 if f[1](x) > 0.
Then h(x)=Heaviside(f[1](x))+2*Heavside(f[2](x)) takes 4 values:

• 0 if f[1](x) < 0 and f[2](x) < 0
• 1 if f[1](x) > 0 and f[2](x) < 0
• 2 if f[1](x) < 0 and f[2](x) > 0
• 3 if f[1](x) > 0 and f[2](x) > 0

Then contourplotting h(x)  with contours=[0, 1, 2, 3] (provided you use a grid dense enough) will display the 4 domains corresponding to each of theconditions above.

I replaced the Heaviside function by a smooth tanh approximation) to ease the computations of the contour levels.
(the smoothing depends on a parameter [set to 1e6] in the tanh.

The generalization of the "Heaviside trick" is

h(x) = add(2^(n-1)*f[n](x), n=1..N)

@Rouben Rostamian  

The values

[E__1 = 0.991324553918355, E__2 = 0.972412189223068, h__1 = 0.999999468863441, nu = 0.473159649082875]

verify eqE.
To get them do

J := (lhs - rhs)(eqE)^2:
opt := Optimization:-Minimize(J, {0 <= nu, nu <= 0.5}, assume = nonnegative)

But I agree that there is probably no solutions:

J := add(`~`[lhs - rhs]([eqA, eqC, eqD, eqE]) ^~ 2);
opt := Optimization:-Minimize(J, {0 <= nu, nu <= 0.5}, assume = nonnegative, iterationlimit = 10000);
       [                      [                          7  
opt := [0.206149604449184010, [E__1 = 3.48734878853157 10 , 

                                                     5         ]]
  E__2 = 13328.4876967435, h__1 = 6.85943362585648 10 , nu = 0.]]


eval([eqA, eqC, eqD, eqE], opt[2]);
              [0.4017000000 = 4.03508624851090, 

                0.1745000000 = -1.93898396374958, 

                0.1517000000 = -1.41034232626533, 

                0.1332000000 = -1.93899197963935]

A simple observation: nu is likely the Poisson coefficient and E1 and E2 are likely Young modulii.`

Thus h__eq (I guess a length?) has the same dimension than h__1.
The relations which define R__C, R__D and R__E do not seem consistent from a dimensional point of view: shouldn't contain h__1^2 instead of h__1

Maybe it could help if you give us the units you use?

More of this the ranges in the fsolve command seem (at least to me) quite weird (if I agree for the nu range, ranges for E are strange).

@sursumCorda 

Thank's a lot.

Edited

Didn't you miss a derivative?
Logistic (ordinary) differential equation writes

diff(u(t), t) = u(t)*(1-u(t))

and Logistic Fractional Equation writes

diff(u(t), t^alpha) = u(t)*(1-u(t))

where  0 < alpha < 1.

I didn't see any derivative in what you wrote

Here is a reference:

https://www.sciencedirect.com/science/article/pii/S0893965921003037

Why don't you ask your question on a Matlab Q&A forum?

@Pepini 

http://math.univ-lyon1.fr/index/homes-www/recrutements/ecoinfo/~borrelli/Articles/PNAS_version_soumise.pdf

The complete construction of the embedding of a flat torus in the 3D Euclidean space is presented here
https://www.emis.de/journals/em/docs/ensaio_matematico/em_24_borrelli-et-al.pdf

Source code avaliable here (hevea.tgz.)
https://hevea-project.fr/ENPageToreCodeSource.html

Good luck if you are courageous enough to translate it in Maple

"I didnt know about Borrelli..."
The Nash-Kuiper embedding theorem(s) has been published in the mid fifties.
This theromem implies the existence of an embedding of a flat torus in the 3D Euclidean space.
Nevertheless, it was not until almost 60 years later that Borrelli, Jabrane, Lazarus and Thibert built a 3D representation of this object.

@Pepini

Look my reply to @Carl Love

@Carl Love

The figure at the end of the question cannot be obtained by rotating the red curve around some axis..

I don't know what is the name of this toroidal figure in English, discovered by the French mathematician Vincent Borrelli (Lyon University, 2012), but it is French dubbed the "Tore plat" (Flat Torus). This "Flat Torus" is an example of smooth fractals.
You can look here
https://www-sop.inria.fr/geometrica/events/wocg14-symmetry_and_periodicity/slides/thibert.pdf

https://www.google.com/search?q=square+flat+torus&tbm=isch&ved=2ahUKEwjQnYnH2q6AAxW9mCcCHcJFBAAQ2-cCegQIABAA&oq=square+flat+torus&gs_lcp=CgNpbWcQAzoFCAAQgAQ6CAgAELEDEIMBOgQIABADOggIABCABBCxAzoHCAAQigUQQzoLCAAQgAQQsQMQgwE6BAgAEB46BwgAEBMQgAQ6CAgAEAUQHhATOggIABAIEB4QE1CRFFiOPWDOQ2gAcAB4AIABcYgB8AiSAQQxNi4ymAEAoAEBqgELZ3dzLXdpei1pbWfAAQE&sclient=img&ei=vErCZNDPHL2xnsEPwosR&bih=808&biw=1729&client=firefox-b-d

This "Flat Torus" is related to the Nash–Kuiper theorem:
https://en.wikipedia.org/wiki/Nash_embedding_theorems

@Rana47 

1. The ode has 3 parameters b, n and k.
    
2. Its HPM counterparty depends on one more parameter wich is the order of the expansion (what I named the ansatz).

Download parameterized_solution.mw

@Rana47 

You write "I have compared the ist and second orders problem but its look not same as reported in paper"
From the paper (equation 33)

• green = zeroth order
• red first order
• blue second order
• 

What I get

• zeroth orders are the same
• first orders are the same given n=5
• second orders differ by a factor n in the paper.
  I noticed this in my mw file sayink that it could be either an error on my side or a typo in the paper.
  Here is the ode for the second order problem from paper ,eq 30.
  Here is this same ode same second in my file is  (remember n=5 and w <--> f, L(w₂)=diff(w₂, x$2))
  it is clear these two odes are identical.
  Given the doundary conditions I use are those of the paper (you can check that easily), and unless Maple is not capable to compute correctly the solution of this problem, I lean towards an error in the paper: the coefficient "n" in the second order term in equation 33 is a typo.
  
  In fact I am even sure there is an error in eq (33)  because after some substitutions the second order term is written this way (eq (34))
  
  As you can see the lultiplicative constant n in the second order term has disappeared.

This should close your first remark.

I dont understand what you say after "Sir, if we suppose ..."
Do you mean that "my" first order problem is not the one given by equations (25)-(26) ???
Does the absence of term k seem like a mistake on my side?
This is normal, this term is balanced by diff(v[0](x), x$2) which is equal to k.

Here is the same file as previously but with notations w instead of f and q instead of p to get closer to to make comparisons easier.
Comparison_made_easier.mw
Read it carefully before saying I did something wrong (which is of course perfectly possible)