@Andiguys 

proposal-2.mw

@dharr 

I spent a lot of time to write my ow, answer to the OP because things that ought to have worked didn't.

Until I noticed that opening the OP worksheet with Maple 2015 reformatted the floats while replacing the decimal point with a comma. Even if some numeric computations remain possible, invoking any plot command generates an error.

Have you ever seen this?
PLOT11_Maple_2015.mw

I add new comments to your previous question

@dharr  @Math-dashti

Nowhere in the paper an excerpt of you provide are mentioned the initial conitions not the integration scheme the authors use.
In page 9/16 a few hints are given: a "Runge-Kutta scheme" (of what order? with wich parameters?) and 4 different initial conditions (none of them corresponding tp those you use). It is imperative to know what initial conditions the authors used.

In the attached file you will find a little Explore(ation) of the system {de1, de2, ic} for theta=1, epsilon in the range 1..3 and
tunable IC:

• 4 sliders (e0, .., e3) enable adjusting with precision the value of epsilon according to the formula
  epsilon : e0 + e1/10 + e2/100 + e3/1000
   
• 3 sliders (u1, .., u3) enable adjusting with precision the value of u(0) according to the formula
  u(0) : u1/10 + u2/100 + u3/1000
   
• 3 sliders (v1, .., v3) enable adjusting with precision the value of v(0) according to the formula
  v(0) : v1/10 + v2/100 + v3/1000

You can play with those sliders and try to find thepatterns depicted in the left column of paper Fig. 2.
chaotic-2_Explore.mw

I haven't been capable to find a pettern ressembling, even from far, to the  paper's Fig. 2 column 1,  row 2.

@dharr @Math-dashti

First point:
Shouldn't de2 write

de2 := diff(v(t),t) = - K[1]*u(t)^3  -  K[2]*u(t) + epsilon*sin(theta*t);

instead of

de2 := diff(v(t),t) = - K[1]*u(t)^3  +  K[2]*u(t) + epsilon*sin(theta*t);

?

Second point:
The differences from between Maple and Matlab could come from the integration method usedn the precisions used... or simply because the figures in the paper you present do not correspond to their captions ot to the eqiuations solved.

The integration method and/or the precision required seem to have only a marginal effect in Maple, leading to infer that these two parameters are likely not the main caise of the Maple-Matlab discrepancy.

Third point:
Based on  @Math-dashti last comment I got the idea ti verify what time step (adaptive of course) rkf45 used.
One can define a new variable to catch this time step and then solve a DAE system instead of a ODE one, but I prefered using the piecewise dsolve output instead (why? I don't know)... and the result was very strange.
As you will see in the next part of the attached file, dsolve with procedure output, even with a highly refined odeplot, gives a quite different result than piecewise dsolve output.

So I wrote a new procedure named newfigures to compare the results its give in the paper Fig.1 row 1 case (only, for matter of time).
Here again the differences between procedure and piecewise outputs are significative (note that piecewise output leads to a very large computing time in comparison to  procedure output).

Burt this induced a very disturbing behaviour of Maple: some remanance effect... look to the the very last part of the attached file and please focus on the red large text 

The worksheet below is based on @dharr's  chaotic-2_mmcdara_4.mw worksheet.
(all plots removed for the file would be to large to upload instead)

chaotic-2_corrected.mw

@dharr 

Thanks for your cooperation and for having finalized what I had started.

I agree with you about the "typo somewhere in the paper": the patterns one may obtain with @Math-dashti primal code are definitely very far from those in the paper.

@dharr 

"If I comment out the range option in dsolve, then using refine=1 in odeplot gives an error:"
Yes, this the same thing in Maple 2015.

"Without the refine option there is insufficient resolution in the plots"
Yes indeed, and this is why I used range+refine.

"Not sure what the solution is."
Do you mean that Maple 2024 cannot reproduce the Maple 2015 results in a simple way?
If so, I suppose you should send an SCR, shouldn't you?

The attached file provides a "smooth" solution without using range+refine: could you you check if it gives better results with Maple 2024 too (TIA)?
chaotic-2_mmcdara_4.mw

@Math-dashti 

You must not simply write Cyan or Gold but "Cyan" or "Gold" (keep the double quotes)

I've just browsed online through the dsolve and odeplot commands for Maple 2024 and haven't found any syntax differences that would explain your empty plots.

I really think you have to wait for someone else to help you.
(sundays are quit calm here, just be patient)

Lastt attempt: replace cyan by "Cyan" and gold by "Gold" when you call figures, indeed Maple 2024 does no longer seem to accept a name for a color and a string seems mandatory.

@Math-dashti 

Can you comment thr line , range=rt[1]..rt[2]  in the dsolve command?

dsol := dsolve(
          {de1,de2,ic}
          , numeric
          , parameters=[epsilon, theta]
      #  , range=rt[1]..rt[2]
          , maxfun=10^5
        );

and tell me if the block labelled "Command 2" still returns a void plot?

@Math-dashti 

Likely a version issue, I use Maple 2015.
Do you get some code error?

It's better for you to wait for someone using Maple 2024 to help you. Meanwhile the code below vould help localizing the issue: can you execute the blocks corresponding to commands 1 to 6 and give me four feedback?

TIA

chaotic-2_mmcdara_2.mw

And a well-designed and instructive serie of examples.
It's not really my line of work and I might never use LAD, but I really enjoyed reading what you did it.
I vote up, of course, twice rather than once

@acer 

sorry for the inconvenience.

@JoyDivisionMan 

This worksheet might onterest you for future needs
Could_be_useful.mw

@acer 

I've tried  this

G := solve({myM=1/2, x>0, x<1},x,explicit)

but I didn't get the single expected solution: why is assuming x>0, x<1 mandatory?
This seems redundant to me.

Just to satisfy mycuriosity: why is the expression produced by myM assuming x>0, x<1 that complex? And why simplifying it provides a factorizes expression whereas I need to use factor(myM2) to get this same form (simplify does nothing)?

Download myM.mw