@tomleslie Yiannis (the OP) is using a 400 x 400 grid and 300 max iterations per point, 88 average iterations per point.

@J4James 

Setting Digits:= 5 gives much faster plots, even at 500..511 t range Preben used. Not every initial condition leads to the cycle (they may lead to singularities, or the at least the numerical method thinks they do), so I chose a random 5. This plot is nearly instantaneous:

randinitial:= combinat:-randcomb(initialset, 5);
Digits:= 5:
DEtools[DEplot](
     [Eq||(1..5)], [x,y,z,v,w](t),
     t= 500..511, maxfun= 0,
     randinitial,
     linecolor= magenta, axes= boxed, scene= [x,y],
     view= [5.2..7.2, 0.2..1.6], scaling= constrained,
     thickness= 0
);

@Preben Alsholm

I noticed that you changed the stepsize from 0.01 to 0.1. This has no effect. The stepsize is only used with classical methods. The default method is rkf45 (same as dsolve), which uses variable step sizes chosen based on error-tolerance settings.

@Chirag I am sorry, but that is far beyond my level of knowledge, and I can't give you any help with it.

@Kitonum 

The algorithm is simply to "walk" along the circle x^2+y^2 = n in a down-the-staircase manner from theta = Pi/4 to theta = 0. If the current point is outside or on the circle, we step down (decrement y); if it is inside the circle, we step forward (increment x). My major improvements were to use isqrt (which is extremely fast) to avoid unit stepping, and to avoid unnecessay recomputation of the squares.

@GPY Certainly. I see how that could be confusing. x2 + 2*x - 1 is the new value of x^2. Here's a proof: Let x be given, let x2= x^2, and let newx= x+1. Then

newx^2 = (x+1)^2 = x^2 + 2*x + 1 = x^2 + 2*x + 2 - 1 = x^2 + 2*(x+1) - 1
           = x2 + 2*newx - 1. Since we've changed x to newx, that becomes x2+2*x-1.

@Chirag 

Yes, of course, if you replace diff(theta(t),t) with theta_prime(t) then theta(t) is completely meaningless to the system. Of course you get 0 as an answer. You need to change the EulerLagrange command to

EulerLagrange(L3, t, theta_prime(t));

Then you will get the answer that you expect.

@J4James 

The complexity of the plot is very affected by your range of t (the independent variable) and by the number of initial conditions that you use. Your original set of initial conditions numbered 180. I decided to work with a randomly selected 5 to 10 of those.

The quality of the plot is greatly increased by setting thickness= 0.

Re-execute this several times to get different random subsets of initial conditions.

randinitial:= combinat:-randcomb(initialset, rand(5..10)());

DEtools[DEplot](
     [Eq1,Eq2,Eq3,Eq4,Eq5], [x(t),y(t),z(t),v(t),w(t)],
     t= 0..140,
     randinitial,
     stepsize= 0.01,
     linecolor= magenta, axes= boxed, scene= [x(t),y(t)],
     view= [5.2..7.2, 0.2..1.6], scaling= constrained,
     thickness= 0
);

@J4James 

If I change the x and y ranges as you have done, I also get an empty plot. I don't know why that is. Perhaps the initial point needs to be within the ranges? If you simply want to crop the plot, use your original ranges and include the option view= [5.2..7.2, 0.2..1.6].

Do you have any Maple code to accompany your Question?

@Jahani_21 

Is this simple enough?

Sin:= proc(x,n)
local k;
     Sum((-1)^k*x^(2*k+1)/(2*k+1)!, k= 0..n-1)
end proc:

@acer 

If you or the OP don't like the summation symbol appearing in gray when using Sum, then you could simply set

`print/Sum`:= ()-> 'sum'(args);

After that, all summation signs should appear the same.

The eigenvalues are purely imaginary and the equilibrium point is a center. This is clearly shown by your phase plot.

@Chia I'm not sure that I understand your question. Are you saying that you want to restrict the solutions to being multiples of 10? Then just use e = 10.

The technique for the second problem is essentially the same as that for the first. The plotting is a little more sophisticated because there are two parameters.

I didn't show the plot for it, but the symmetrically corresponding region below the horizontal axis is the region for which (0,0,0) is a relative maximum.

Download 2nd_derivative_test_2.mw