@mmcdara The worksheet opens with the comma's but executing it changes them to ".", but the plot error is the same owing to complex values. 

@Math-dashti I don't understand exactly what you want - please specify in the notation of your worksheet exactly what you want. I think S__1 and S__2 are some functions of p1,p2,p3,k? If so I don't know what  functions those are and what you want to do with them?

Or if you want to reproduce the paper results, please give a worksheet with that notation.

@Andiguys As I said, the view option is the syntax to magnify it. but it isn't helpful in this case.

@Math-dashti S__1 and S__2 are not mentioned in the picture you supplied so I do not know what you are asking.

Yes, you cannot plot complex values on a real plot. Don't you want to reject complex equilibria? Aren't only real u(t), v(t) of interest?

@mmcdara  I tried to reproduce the behaviour in a minimal example for an scr, but things are different in 2025 from 2024, so I am giving up on this.

@Math-dashti For a chaotic system we may expect the solution to be very sensitive to the type of solver and step size algorithm so it is not surprising that Matlab and Maple give somewhat different results.

Edit: I see you used different parameters for the two cases. But I can't reproduce Fig 1a even approximately, so I think there is a typo somewhere in the paper.

@mmcdara Your chaotic-2_mmcdara_4.mw works in 2024:

chaotic-2_mmcdara_4_2024.mw

In the earlier behaviour, there was a complicated interaction where range+refine worked the solution up to a point, then suggested increasing maxfun to go beyond that point, but then gave a warning and didn't work if maxfun was too high. But without parameters it worked with those same settings (see other reply). Range and refine work with other ODEs, e.g. the Van der Pol one in the odeplot help page, so it is something about accuracy, which you solved with relerr.

@mmcdara @Math-dashti ​​​​​​​Here is a working Maple 2024 version of @mmcdara's code. The workaround was to remove the parameters functionality of dsolve, and just do a new dsolve each time.

chaotic-2_mmcdara_4.mw

@mmcdara (in Maple 2024) If I comment out the range option in dsolve, then using refine=1 in odeplot gives an error:
"Error, (in plots/odeplot) 'refine' option can only be used for the procedure output of the 'rkf45', 'ck45', 'forwardeuler', 'backwardeuler', or 'rosenbrock' methods used with the 'range' argument"

Without the refine option there is insufficient resolution in the plots. Not sure what the solution is.
The color specifications are not a problem. Here without range or refine:

chaotic-2_mmcdara_3.mw

@Math-dashti You start by setting the coefficients in E8 to zero. But if you did that with E9 you would find from the coefficient of U^3 that k[2]=0. Is that allowed or is k[2] a given? I do not think it is correct to set the coefficients equal, e.g., in W1, becase both E8 and E9 are equal to zero, so you could multiply them by anything before you set them equal and then would get a different answer. The elimination method I used amounts to multiplying them by something that makes the coefficient of y'' the same, and then comparing the coefficients - if I do this manually I get the same result (see attached). Perhaps I should use the other condition to say, eliminate zeta[1] or zeta[2]? but they both appear in the desired result.

Download C2.mw

@Math-dashti It seems DEplot isn't set up for this. I'm not sure about fieldplot. Of course you could convert back to cartesian coordinates and use DEplot. odeplot works but is tedious.

Download polarphase.mw

@Math-dashti The output from PDEtools:-ConservedCurrents({de1,de2}, [u(t),v(t)]); correctly says any function of it will do. I stripped off the arbitrary function and just multiplied by that constant just to agree with @Rouben Rostamian's result. If something is conserved, then any function of it is also conserved.

@Math-dashti FYI, Maple can calculate the conserved quantity from the odes by

PDEtools:-ConservedCurrents({de1,de2}, [u(t),v(t)]);
-(1/2)*op(1, rhs(op(%)));

(The -1/2 just to convert to @Rouben Rostamian's form.)

@Rouben Rostamian  I'm fascinated by symbolic solutions, which was my main interest. Otherwise, I'm not convinced it is that useful. There have been many requests about it on Mapleprimes, so it was more a coding challenge.

@Alfred_F Here's one way. You always have to play around with these things.

Download test.mw