restart:

kernelopts(version);

Procedure

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*alpha/(4*k__1*q^2/d^2-alpha__3*xi/eta__1 - chi*H^2),
                       n:= 50,
                       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,soln3, Imagroot1, Imagroot2, AreTherePurelyImaginaryRoots;

# Assign g for q and plot g, Set q as a complex and Compute the Special Asymptotes

  g:= q-(1-alpha)*tan(q)+  c*tan(q) + w*q:
  f:= subs(q = x+I*y, g):
  b1:= evalc(Re(f)) = 0:
  b2:= y-(1-alpha)*tanh(y) -(alpha__3*xi*alpha/(eta__1*(4*k__1*y^2/d^2+alpha__3*xi/eta__1)))*tanh(y) = 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 complex roots

  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;

AreTherePurelyImaginaryRoots:= type(true, 'truefalse');
  try
# Compute the rest of the Roots
  soln3:= simplify~(fsolve({ b2}, { y = 0 .. infinity}),zero);  
  Imagroot1:=subs(soln3, I*y);
  Imagroot2:= -Imagroot1;
  catch:
   AreTherePurelyImaginaryRoots:= type(FAIL, 'truefalse');
  end try;

## odd and all asymptotes
  OddAsymptotes:= Vector[row]([seq(evalf((1/2)*(2*j+1)*Pi), j = 0 .. n)]);
  AllAsymptotes:= sort(Vector[row]([OddAsymptotes, qstar]));
       
 if (xi > 0) then
  if AreThereComplexRoots and not AreTherePurelyImaginaryRoots
  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 AreThereComplexRoots and AreTherePurelyImaginaryRoots
  then gg:= [qcomplex1, qcomplex2, Imagroot2, 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 and AreTherePurelyImaginaryRoots
  then gg:= [Imagroot2, 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 and not AreTherePurelyImaginaryRoots
  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:
else
if AreThereComplexRoots and not AreTherePurelyImaginaryRoots
  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 AreThereComplexRoots and AreTherePurelyImaginaryRoots
  then gg:= [qcomplex1, qcomplex2, Imagroot2, 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 and AreTherePurelyImaginaryRoots
  then gg:= [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 and not AreTherePurelyImaginaryRoots
  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:
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] := simplify(evalf(gamma__1*alpha/(4*k__1*qq[i]^2/d^2 - alpha__3*xi/eta__1- chi*H^2)),':-zero'):
   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):
  r := LinearSolve(A,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):

## compute minimum value of tau
  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,qq;
     end proc:

Run the calculation for supplied value of 'xi

ans:=[doCalc(0.1, 2)]:
evalf(ans[9]);
evalf(ans[10]);

Contour plot for pmin (minimum  of tau) for different xi and H values