13202 Reputation

12 years, 350 days

Recovering data...

Assuming you have the automatic back-up on, ie the menu entry

Tools->Options->General->Auto save entry

checkbox is "checked",  then Maple will be backing up your work at the interval specified. These backups can be recovered in a susequent Maple session using

File->Restore Backup

which will produce a pop-up window where you can browse for the backup you want. (Rather amusingly, in this pop-up, one has to uses the horizontal scroll bar in order to scroll vertically!). Alternatively, you can use the

File->Open

and navigate  to the backup directory, which in a standard Windows installlation is located at

A few warnings

1. If the checkbox Tools->Options->General->Keep files is set (the default) then you will keep the backup files even after using a "Save" or "Save as" command. If this is unset, then the "Save" or "Save as" command, performs the "save" operation and deletes the backup.
2. It is generally good practice when you open a new Maple worksheet to save it immediately with "sensible" name. Otherwise, (assuming Keep files checkbox is set), your backup directory will contain a lot of files called something like "untitled3214_MAS.bak" - which isn't exactly helpful. You can, of course, generally figure out which one you want from the datestamp/timestamp.
3. The default interval (which you can change) between backups is 3minutes. If you have a BIG worksheet, then you may notice that Maple becomes "unresponsive" for a few seconds every 3mins. This is unavoidable. I'm pretty sure that when Maple is performing a backup, the annunciator in the bottom left corner of the Maple GUI, says "saving" - so at least you can check that Maple hasn't died on you! Increasing the back-up interval is a trade-off between how much work you risk losing, and how annoying you find the periodic non-response.

OOps...

Just noticed that you want this evaluated for m=1..10 - soI need  a trivial additional execution group to my previous post. See the attached

 >
 (1)
 >
 (2)
 >

A suggestion...

This problem is (somewhat?) analogous to a Sudoku number puzzle, in that it is grid-based and there are certain constraints on what has to occur in each row and column.

There is a wonderful piece of code here (https://www.maplesoft.com/Applications/Detail.aspx?id=154483), which uses Maple's Logic:-Satisfy() command to solve "the world's hardest sudoku".

If one looks up the help for the Logic:-Satisfy() command one might get the impression that it is applicable to only rather "simple"  logical problems. However in the Application mentioned above, the author generates 17598 logical constraints to supply to the Logic:-Satisfy() command, which then solves the Sudoku puzzle in ~0.5secs.

I'm pretty sure that a similar approach, with the appropriate constraints, could be used for the OP's problem.

Do I fancy writing the necessary constraints? Not really!  I think this could turn into a rather lengthy programming exercise

Now concatention makes sense...

but only for the labels of each point. See the attached. (Again, the plot renders better in the worksheet than it does on this site)

 > restart:   with(plots):   theta := (2*Pi)/17:   p1:= pointplot        ( [ seq            ( [cos(k*theta), sin(k*theta)],               k=1..16            )          ],          symbol = solidcircle,          color = red,          symbolsize = 10        ):   p2:= textplot        ( [ seq            ( [ cos(k*theta), sin(k*theta), cat("M", k)],              k=1..16            )          ],          align=["above", "right"]        ):   display([p1,p2]);
 >

Use the big green up-arrow...

In the Mapleprimes toolbar to upload the actual worksheet you are running which illustrates the problem.

Well...

the first thing you can do, is to utilise the amzing featued of this site where you can upload the worksheet which is causing you problems. This is done using the big green up-arrow in the Mapleprimes toolbar

When you cut/paste from a Maple worksheet and expect people to trawl throught it deleting all the irrelevant (output) crap - you are on a loser - opload something executable which illustrates your problem

Well...

These mean two completely different things and I'm not sure what you are after.

p(Delta*delta)  is the function p() with the argument (delta*delta)

p*(Delta*delta)is just a simple product of three terms, and the parentheses are completely superfluous

Only you know whether you have a function definition or a simple product.

Since have not...

uploaded a worksheet (use the big green arrow in the Mapleprimes toolbar), it is difficult to work out exactly what your problem is. Maybe the attached will help?

 > restart:   with(PDEtools): # # The "classic" wave equation #   we1:= diff(u(x,t), t$2)= c^2*diff(u(x,t), x$2); # # The "usual" variable transformation for the # above #   we2:= dchange         ( solve           ( [ eta=x-c*t, mu=x+c*t],             {x,t}           ),           we1,           {eta, mu},           params={c}         ):   we3:=simplify( rhs(we2)-lhs(we2)=0); # # Now solve we3 #   pdsolve(we3);
 (1)
 > # # Note that one can get the same answer more or less directly #   pdsolve(we1);
 (2)
 >

There are four...

main plot commands which provide "contours". These are

contourplot()

contoutplot3d()

listcontplot()

listcontplot3d()

The first two of these produce contours for continuous expressions: the second two produce contours for a "grid" of discrete data. So far as I am aware, you only have data at discrete points, so you are restricted to the last two entries above. You may be able to restructure your code (slightly?) to use either of the first two, if you define a procedure whihc returns a meaningful answer for a pair of supplied variables whose ranges you provide. You (probably?) won't be able to specify the exact points which will be used in this case, since these will be determined by Maple's adaptive plotting routines.

See the attached "toy" example for usage, where I produce plots for the "same" expression using all four commands. (I haven't bothered to make these look as similar as possible by adjusting color schemes, number of contours, axes positions, etc, etc)

 > restart;   with(plots):   f:=(x,y)-> sin(x/5*y/5); # # Produce a standard 2D contour plot of the # above expression #   contourplot( f, -15..15, -15..15, filledregions=true); # # Produce a standard 3D contour plot of the # above expression #   contourplot3d( f, -15..15, -15..15, filledregions=true); # # Sample the above expression on a 31 by 31 grid to a listlist # structure and use listcontplot to generate a 2D contourplot # of the result #   listcontplot([seq([seq( f(x,y), x=-15..15)], y=-15..15)], filledregions=true); # # Sample the above expression on a 31 by 31 grid to a listlist # structure and use listcontplot3d to generate a 3D contourplot # of the result #   listcontplot3d([seq([seq( f(x,y), x=-15..15)], y=-15..15)], filledregions=true);
 >

You...

@ijuptilk

generate 3 lists each of which has 21 entries so plot axes will be 1..3 nd 1..21. If this isn't clear enough consider the attached "toy" example from the listcontplot help pags

 > restart;
 > with(plots):
 > #   # there are 31 lists with 31 entries, so both axes   # will run from 1 to 31.   #     listcontplot([seq([seq( sin(x/5*y/5), x=-15..15)], y=-15..15)]);
 > #   # there are 11 lists with 31 entries, so one axis   # will be 1..11 and the other will be 1..31   #     listcontplot([seq([seq( sin(x/5*y/5), x=-15..15)], y=-5..5)]);
 > #   # there are 11 lists with 7 entries, so one axis   # will be 1..11 and the other will be 1..71   #     listcontplot([seq([seq( sin(x/5*y/5), x=-3..3)], y=-5..5)]);
 >

Rarely...

have I seen code with so many basic logical, mathematical and syntactical errors.

In the attached I attempted to list these and suggest corrections based on what I guessed to be your intention, but eventually I had to give up.

So for what it is worth

 > restart: # # OP's attempt at code # #     with(PDEtools):L:=2: #     PDEtools[declare](theta(r, z), prime = r, prime = z); #     PDEtools[declare](sigma(r, z), prime = r, prime = z); # # Problems: # # 1. Why load the PDEtools package and the use the longform #    calling sequence PDEtools[declare]() ? # 2. 'prime' notation can only be used for functions of one #    variable. If I were to write g'(x,y), with respect to #    which variable has the differentiation been performed? # # Since OP obviously does not understand clever notation, or # how to use it, probably safest not to bother with it, in # which case the above attempt at code becomes #   L:=2:
 > # # OP's next attempt at code # #      theta(r, z):=sum((p^i)*theta[i](r, z),i=0..L):   #      sigma(r, z):=sum((p^i)*sigma[i](r, z),i=0..L); # # Problems: # # 1. *Never* use an indexed variable, and an unindexed variable #    with the saem name. For exampl theta[i](r,z) in the above #    defines the i-th entry in a table named theta: Is the theta #    on the left-hand side of this assignment the same table? Or #    maybe something "different. This confuses me, and it confuses #    Maple # 2. Don't use sum() when you should be using add() # # Fix the above attempt at code - note the capitalization on # the left-hand side of the assignments to avoid the name clash # described in (1) aabove #   Theta(r, z):= add((p^i)*theta[i](r, z), i=0..L);   Sigma(r, z):= add((p^i)*sigma[i](r, z), i=0..L);
 (1)
 > # # OP's next attempt at code # # HO1:=(1-p)*(diff(theta(r, z),r,r))+p*(((1/r)*(diff(theta(r, z),r))+Nb*((diff(theta(r, z),r))*diff(sigma(r, z),r))+Nt*(diff(theta(r, z),r)^2))): # expand(%): # collect(%,p): # HO2:=(1-p)*(diff(sigma(r, z),r,r))+p*(((1/r)*(diff(sigma(r, z),r))+(Nt/Nb)*((diff(delta(r, z),r,r))+(1/r)*diff(delta(r, z),r)))): # expand(%): # collect(%,p): # # # Problems. # # 1. All instances of 'theta(r,z) and 'sigma(r,z) in the above will #    have to be replaced with 'Theta(r,z) and 'Sigma(r,z)' for reasons #    described in the comments of the preceding execution group # 2. The expand() and collect() function in the above are pointless. #    Their output is not assigned to anything and so cannot be #    subsequently called by anything. # 3. OP has two equations, but even after collecting coefficients of #    powers of 'p' later ther will be *three* dependent variables. #    This is not going to end well, but leave it for now! #    Where did delta(r,z) come from? #   HO1:=(1-p)*(diff(Theta(r, z),r,r))+p*(((1/r)*(diff(Theta(r, z),r))+Nb*((diff(Theta(r, z),r))*diff(Sigma(r, z),r))+Nt*(diff(Theta(r, z),r)^2)));   HO2:=(1-p)*(diff(Sigma(r, z),r,r))+p*(((1/r)*(diff(Sigma(r, z),r))+(Nt/Nb)*((diff(delta(r, z),r,r))+(1/r)*diff(delta(r, z),r))));
 (2)

 > # # OP's next code attempt # #     for i from 0 to L+1 do equa[1][i]:=coeff(HO1,p,i)=0 end  do: #     for i from 0 to L+1 do equa[2][i]:=coeff(HO2,p,i)=0 end  do: # # Questions # # 1. Why use two 'for' loops when one will suffice? # 2. Why use indexed variable on the left-hand-side #    of these assignments when there is no need. This #    just makes "notation" even more confusing #   for i from 0 to L+1 do       equa1[i]:=coeff(HO1,p,i)=0;       equa2[i]:=coeff(HO2,p,i)=0   end  do;
 (3)
 > # # OP's next code attempt, even more hilarious # #    con[1][0]:=theta(r,z)[0](h(z),0)=0,(D(theta(r,z)[0]))(0,0)=0: #    con[2][0]:=sigma(r,z)[0](h(z),0)=0,(D(sigma(r,z)[0]))(0,0)=0: #    for j from 1 to L do: #         con[1][j]:=theta(r,z)[j](h(z),0)=0,(D(theta(r,z)[j]))(0)=0: #         con[2][j]:=sigma(r,z)[j](h(z),0)=0,(D(sigma(r,z)[j]))(0)=0: #    end do; # # Problems # # 1. There are no functions theta(r,z)[0] or sigma(r,z)[0], these #    should presumably be theta[0](r,z) and sigma[0](r,z) # 2. The correction in (1) above shows just how preposterous #    the original definition is because then we would have #    theta[0](r,z)(h(z), 0) - which is a "sensible" function #    isn't it!!! # 3. The derivative conditions above have the same inept function #    definition. Furthermore the 'D' operator has no indication #    of which independent variable differetniation is with respect #    to. In the corrected code below I have assume that it is the #    first variable # 4. Why is the case j=0 treated separately from the others? There #    seems to be no reason! Just run the loop from 0 to L. # 5. Why the unnecessary indexing on the left hand side of the #    assignements - again just to make the notation as confusing #    as possible? #   for j from 0 to L do       con1[j]:= theta[j](h(z),0)=0, D[1](theta[j])(0,0)=0:       con2[j]:= sigma[j](h(z),0)=0, D[1](sigma[j])(0,0)=0:   end do;
 (4)
 > # # This is when I lost the will to live # #       for i from 0 to L do: #           dsolve({equa[1][i],con[1][i]},delta(r,z)[i](r)): #           delta(r,z)[i](r):=rhs(%); #           delta(r,z)[i](r):=evalf(%); #           dsolve({equa[2][i],con[2][i]},sigma(r,z)[i](r)): #           sigma(r,z)[i](r):=rhs(%); #           sigma(r,z)[i](r):=evalf(%); #       end do; # # There are so many errors (logical and syntactical) in the # above that I don't even know where to start # # 1. What does this delta(r,z)[i](r) mean? A function #    theta() which takes two variables (r,z) which returns an #    indexable quantity 'delta(r,z)[i]' which in trurn is a function #    function accepting the single variable r ???? # 2. Ditto for the gobbledegook sigma(r,z)[i](r). # 3. So the OP apparently wants to solve differential equations #    for functions whihc do not exist in the differetnial equations # # More out of idle curiosity rather than anything else. what is # the system to be solved first time the loop executes (ie i=0) #   eval( [equa1[0], con1[0]], i=0); # # And dsolve() will Error on this because of the invalid boundary # condition involving the unknown function h(z) #   dsolve(%); # # Check the other equation #   eval( [equa2[0], con2[0]], i=0); # # dsolve() will error on this for the same reason as above
 (5)
 >

your latest post very carefully, I now have absolutely no idea what you are trying to achieve.

The secret to getting a sensible answer is to ask a sensible question

Yup...

I get the same error as you do in Maple 2021. Looks like a bug that has been fixed

You state...

why do I have repeated values for j and u? And how can solve this problem without using Makeunique function.

Thhis is because the loop system I set up treats 'u' and 'j' as independent variables, so in the worksheet I provided initially, the command

ans1:= CodeTools:-Usage(Array([seq(seq( [j, doCalc(j,u)], j=-2..0, 2), u =1..334,333)])):

sets u=1, then does the calculation for j=-2 and 0, the sets u=334 amd does the calculation for j=-2 and 0. It is still not clear to me what you actually want to achieve, unless 'j' is not an independent variable, but is in fact a simple function 'u'. One possibility is shown in the attached where as u varies from 1..334, j varies from -2..0, according to j=-2.0+(u-1)/334.

(In the attached, purely to save execution time I do this calculation with 'u' varying in steps of 33. If you want to do this caclulation ffor all u-values from 1..334, then just change 'ustep' to 1

#######################################################

I also want to do contourplot for u vs j for complex values of p_min.

So far as I can tell all returned "values" of p_min are expressions containing the variables 'x' and 't'. Maybe some of these expression will be complex for somevlaues of 'x' and 't', and maybe not. Since I have absoutely no idea what values of 'x' and 't' you paln on using, I haven't really considered this question further

 > restart:   doCalc:= proc( xi , u)  #u is the \bar(H): normalize magnetic field magnitude,                           # where H = bar(H)*H__a                  # Import Packages                  uses ArrayTools, Student:-Calculus1, LinearAlgebra,                       ListTools, RootFinding, plots, ListTools:                  local gamma__1:= .1093,                        alpha__3:= -0.1104e-2,                        k__1:= 6*10^(-12),                        d:= 0.2e-3,                        theta0:= 0.1e-3,                        eta__1:= 0.240e-1, chi:= 1.219*10^(-6),                        alpha:= 1-alpha__3^2/(gamma__1*eta__1),                        theta_init:= theta0*sin(Pi*z/d),                        H__a:= Pi*sqrt(k__1/chi)/d,                        H := u*H__a,                             c:=alpha__3*xi*alpha/(eta__1*(4*k__1*q^2/d^2-alpha__3*xi/eta__1 - chi*H^2)),                        w := chi*H^2*eta__1*alpha/(4*k__1*q^2/d^2-alpha__3*xi/eta__1 - chi*H^2),                        n:= 20,                        g, f, b1, b2, qstar, OddAsymptotes, ModifiedOddAsym,                        qstarTemporary, indexOfqstar2, qstar2, AreThereComplexRoots,                        soln1, soln2, qcomplex1, qcomplex2, gg, qq, m, pp, j, i,                        AllAsymptotes, p, Efun, b, aa, F, A, B, Ainv, r, theta_sol, v, Vfun, v_sol,minp,nstar; # Assign g for q and plot g, Set q as a complex and Compute the Special Asymptotes                  g:= q-(1-alpha)*tan(q)- (w*q + c)*tan(q):                  f:= subs(q = x+I*y, g):                  b1:= evalc(Re(f)) = 0:                  b2:= evalc(Im(f)) = 0:                  qstar:= (fsolve(1/c = 0, q = 0 .. infinity)):                  OddAsymptotes:= Vector[row]([seq(evalf(1/2*(2*j + 1)*Pi), j = 0 .. n)]); # Compute Odd asymptote                  ModifiedOddAsym:= abs(-~(OddAsymptotes, qstar));                  qstarTemporary:= min(ModifiedOddAsym);                  indexOfqstar2:= SearchAll(qstarTemporary, ModifiedOddAsym);                  qstar2:= OddAsymptotes(indexOfqstar2); # Compute Odd asymptote                  AreThereComplexRoots:= type(true, 'truefalse');                  try                       soln1:= fsolve({b1, b2}, {x = min(qstar2, qstar) .. max(qstar2, qstar), y = 0 .. infinity});                       soln2:= fsolve({b1, b2}, {x = min(qstar2, qstar) .. max(qstar2, qstar), y = -infinity .. 0});                       qcomplex1:= subs(soln1, x+I*y);                       qcomplex2:= subs(soln2, x+I*y);                  catch:                        AreThereComplexRoots:= type(FAIL, 'truefalse');                  end try; # Compute the rest of the Roots                  OddAsymptotes:= Vector[row]([seq(evalf((1/2)*(2*j+1)*Pi), j = 0 .. n)]);                  AllAsymptotes:= sort(Vector[row]([OddAsymptotes, qstar]));                  if   AreThereComplexRoots                  then gg:= [ qcomplex1,                              qcomplex2,                              op(Roots(g, q = 0.1e-3 .. AllAsymptotes[1])),                              seq(op(Roots(g, q = AllAsymptotes[i] .. AllAsymptotes[i+1])), i = 1 .. n)                           ];                  elif not AreThereComplexRoots                  then gg:= [ op(Roots(g, q = 0.1e-3 .. AllAsymptotes[1])),                              seq(op(Roots(g, q = AllAsymptotes[i] .. AllAsymptotes[i+1])), i = 1 .. n)                            ];                  end if: # Remove the repeated roots if any & Redefine n                  qq:= MakeUnique(gg):                  m:= numelems(qq): ## Compute the time constants                  for i to m do                  p[i] := gamma__1*alpha/(4*k__1*qq[i]^2/d^2 - alpha__3*xi/eta__1- chi*H^2):                  end do: ## Compute θ_n from initial conditions                 for i to m do                 Efun[i] := cos(qq[i]-2*qq[i]*z/d)-cos(qq[i]):                 end do: ## Compute integral coefficients for off-diagonal elements θ_n matrix                 printlevel := 2:                 for i to m do                     for j from i+1 to m do                         b[i, j] := int(Efun[i]*Efun[j], z = 0 .. d):                         b[j, i] := b[i, j]:                         aa[i, j] := b[i, j]:                         aa[j, i] := b[j, i]:                     end do :                 end do: ## Calculate integral coefficients for diagonal elements in theta_n matrix                 for i to m do                 aa[i, i] := int(Efun[i]^2, z = 0 .. d):                 end do: ## Calculate integrals for RHS vector                for i to m do                F[i] := int(theta_init*Efun[i], z = 0 .. d):                end do: ## Create matrix A and RHS vector B                A := Matrix([seq([seq(aa[i, j], i = 1 .. m)], j = 1 .. m)]):                B := Vector([seq(F[i], i = 1 .. m)]): ## Calculate inverse of A and solve A*x=B               Ainv := 1/A:               r := MatrixVectorMultiply(Ainv, B): ## Define Theta(z,t)              theta_sol := add(r[i]*Efun[i]*exp(-t/p[i]), i = 1 .. m): ## Compute v_n for n times constant              for i to m do              v[i] := (-2*k__1*alpha__3*qq[i])*(1/(d^2*eta__1*alpha*gamma__1))+ alpha__3^2*xi/(2*(eta__1)^2*qq[i]*alpha*gamma__1)+xi/(2*eta__1*qq[i]) + alpha__3*chi*H^2/(2*eta__1*qq[i]*gamma__1*alpha):              end do: ## Compute v(z,t) from initial conditions              for i to m do              Vfun[i] := d*sin(qq[i]-2*qq[i]*z/d)+(2*z-d)*sin(qq[i]):              end do: ## Define v(z,t)              v_sol := add(v[i]*r[i]*Vfun[i]*exp(-t/p[i]), i = 1 .. m): ##              minp:=min( seq( Re(p[i]), i=1..m) ):              member(min(seq( Re(p[i]), i=1..m)),[seq( Re(p[i]), i=1..m)],'nstar'): ## Return all the plots                  return theta_init, theta_sol, v_sol, minp, eval(p), nstar, p[nstar], g, H, H__a;                  end proc:
 > # Run the calculation for supplied value of 'xi'
 > # Put the returned quantities in a simple list
 > ans:=[doCalc(-0.06, 0.001)]: ans1:=[doCalc(0.06, 0.001)]: evalf(ans[9]); evalf(ans[10]);
 (1)
 > ##
 > with(plots): d:= 0.2e-3: ## director Plot negative activity display( plot(ans[1], z=0..d, color = green),          plot([seq(eval(ans[2], t = j), j in [seq(.1*j, j = 1 .. 5)])], z=0..d, legendstyle = [font                 = ["ROMAN", 12],location = right], labeldirections = ["horizontal", "vertical"]),                     axis[1]=[tickmarks=[0="0", 0.00005="50", 0.00010="100", 0.00015="150", 0.00020=                         "200"]], axis[2]=[tickmarks=[0="0", 0.00002="2e-05", 0.00004="4e-05", 0.00006=                     "6e-05", 0.00008= "8e-05", 0.0001="1e-04"]], axes=boxed, labels =[typeset(z," ()                     "), typeset(theta(z, t)," (degrees)")], labelfont = ["TimesNewRoman", 14],                             axesfont = [14, 14]); ## director angle Plot postive activity display( plot(ans1[1], z=0..d, color = green, legend = theta[0]),          plot([seq(eval(ans1[2], t = j), j in [seq(.1*j, j = 1 .. 5)])], z=0..d, legend = [t[1] =                  .2, t[2] = .4, t[3] = .6, t[4] = .8, t[5] = 1.0], legendstyle = [font = ["ROMAN", 12],                location = right], labeldirections = ["horizontal", "vertical"]), axis[1]=[tickmarks=                [0="0", 0.00005="50", 0.00010="100", 0.00015="150", 0.00020="200"]], axis[2]=                        [tickmarks=[0="0", 0.00002="2e-05", 0.00004="4e-05", 0.00006="6e-05", 0.00008=                       "8e-05", 0.0001="1e-04"]], axes=boxed, labels =[typeset(z," ()"), typeset(theta(z, t)," (degrees)")], labelfont =["TimesNewRoman", 14], axesfont = [14, 14]); ## v Plot negative activity plt2:= plot([seq(eval(ans[3], t = j), j in [seq(.2*j, j = 1 .. 5)])], z = 0 .. d, legendstyle = [font = ["ROMAN", 12],                location = right], labeldirections = ["horizontal", "vertical"]): display( plt2, axis[1]=[tickmarks=[0="0", 0.00005="50", 0.00010="100", 0.00015="150", 0.00020="200"]          ], axes=boxed, labels =[typeset(z," ()"), typeset(v(z, t)," (m/s)")], labelfont = ["TimesNewRoman", 14], axesfont = [14, 14]);   ## v Plot postive activity plt21:= plot([seq(eval(ans1[3], t = j), j in [seq(.2*j, j = 1 .. 5)])], z = 0 .. d, legend = [t[1] =                .2, t[2] = .4, t[3] = .6, t[4] = .8, t[5] = 1.0], legendstyle = [font = ["ROMAN", 12],                location = right], labeldirections = ["horizontal", "vertical"]): display( plt21,axis[1]=[tickmarks=[0="0", 0.00005="50", 0.00010="100", 0.00015="150", 0.00020="200"]          ], axes=boxed, labels =[typeset(z," ()"), typeset(v(z, t)," (m/s)")], labelfont = ["TimesNewRoman", 14], axesfont = [14, 14]);
 > # q plot
 > #
 > plot(ans[8], q=0..5*Pi,view=[DEFAULT,-10..10], labels =[q, r(q)], labelfont = ["TimesNewRoman", 14], labeldirections = ["horizontal", "vertical"] );
 > # p values
 > evalf(ans[5]);
 (2)
 > # Director angle plot for at d/2 as a funtion of time
 > theta_f := subs(z=d/2, Re(ans[2])): plot(theta_f, t=0..1000, view=[DEFAULT, DEFAULT], labels =[typeset(t," (seconds)"), typeset(theta('d/2', t),"(rad)")], labelfont = ["TimesNewRoman", 14], axesfont = [14, 14], labeldirections = ["horizontal", "vertical"], axes=boxed); theta_f := subs(z=d/2, ans1[2]): plot(theta_f, t=0..1, view=[DEFAULT,DEFAULT], labels =[typeset(t," (seconds)"), typeset(theta('d/2', t),"(rad)")], labelfont = ["TimesNewRoman", 14], axesfont = [14, 14], labeldirections = ["horizontal", "vertical"], axes=boxed);
 > #Convert and store xi and doCalc into an array   #ans1:= Array([seq( [j, doCalc(j)], j=-2..0, 0.006)]):  ustep:=33:  ans1:= CodeTools:-Usage(Array([seq([-2.0+(u-1)/334, u, doCalc(-2.0+(u-1)/334, u)], u =1..334, ustep)])):
 memory used=4.32GiB, alloc change=32.00MiB, cpu time=47.36s, real time=43.19s, gc time=5.96s
 > #
 > #
 > p_min:= convert(ans1[..,5],list): indets~(p_min, name); j_value := convert(ans1[..,1],list); u_value := convert(ans1[..,2],list); # plot( convert(ans1[..,1],list), p_min/30) eval(p_min[1], [z=1, t=1]);
 (3)
 > ## please why do I have repeated values for j and u? And how can solve this problem without using Makeunique function.
 > ## I also want to do contourplot for u vs j for complex values of p_min.

You further state

Cannot reproduce...

This appears to work for me, see the attached. It may be a version issue: I'm running Maple 2022, wh8hc Maple version are you running?

 > restart:   interface(version);   with(GroupTheory):   CompositionSeries(SmallGroup(336,209));
 (1)
 >