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]);

  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);

#
# Note that one can get the same answer more or less directly
#
  pdsolve(we1);

  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);

  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)]);

  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);

#
# 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))));

#
# 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;

#
# 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;

#
# 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

Error, (in dsolve) unexpected occurrence of the variables {z} in the 1st operand of theta[0](h(z),0) in the given initial conditions

  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]);

##

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]);

  

Error, missing operator or `;`

# 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]);

 # 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]);

## 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.

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

  restart:
  eqn1 := diff(f(eta), `$`(eta, 3))+f(eta)*(diff(f(eta), `$`(eta, 2)))-(diff(f(eta), eta))^2-lambda*(diff(f(eta), eta)) = 0:
  eqn2 := diff(theta(eta), `$`(eta, 2))+f(eta)*(diff(theta(eta), eta))*Pr = 0:
  Bcs := (D(f))(0) = 1, f(0) = 0, (D(f))(infinity) = 0, theta(0) = 1, theta(infinity) = 0:
  params:=[lambda = .5, Pr = 6.3, infinity=10]:
  sol:=dsolve( eval([eqn1, eqn2, Bcs], params), numeric):
  plots:-odeplot( sol,
                  [ [ eta, f(eta) ],
                    [ eta, theta(eta) ]
                  ],
                  eta=0..10,
                  color=[red, blue]
                );