restart:

doCalc:= proc(xi, u)  #u is the \bar(H): normalize magnetic field magnitude,
                          # where H = bar(H)*H__a
                 option remember;
                 # 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;

                 if not [xi,u]::[numeric,numeric] then
                    return 'procname'(args);
                 end if;

# 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(evalhf(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(evalhf((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[n], numeric))
                            ];
                 elif not AreThereComplexRoots
                 then gg:= [ op(Roots(g, q = 0.1e-3 .. AllAsymptotes[n], numeric))
                           ];
                 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] := evalf(gamma__1*alpha/(4*k__1*qq[i]^2/d^2 - alpha__3*xi/eta__1- chi*H^2)):
if not p[i]::complex(numeric) then print(p[i], [xi,u], qq[i]); end if;
                 end do:

                Digits := 15;
## 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] := evalf(Int(Efun[i]*Efun[j], z = 0 .. d)):
                        if not b[i, j]::complex(numeric) then
                            b[i, j] := evalf[15](Int(Efun[i]*Efun[j], z = 0 .. d,
                                                 digits=15, method=_Dexp, epsilon=1e-12)):
                            if not b[i, j]::complex(numeric) then
                               print("A",[Efun[i]*Efun[j], z = 0 .. d]);
                            end if;
                        end if;
                        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] := evalf(Int(Efun[i]^2, z = 0 .. d)):
                   if not aa[i, i]::complex(numeric) then
                      aa[i, i] := evalf[15](Int(Efun[i]^2, z = 0 .. d,
                                                digits=15, epsilon=1e-12));
                      if not aa[i, i]::complex(numeric) then
                         print("B",[Efun[i]^2, z = 0 .. d]);
                      end if;
                   end if;
                end do:

## Calculate integrals for RHS vector

               for i to m do
               F[i] := evalf(Int(theta_init*Efun[i], z = 0 .. d));
               if not F[i]::complex(numeric) then
                  F[i] := evalf[15](Int(theta_init*Efun[i], z = 0 .. d,
                                     digits=15, method=_Dexp, epsilon=1e-12));
                  if not F[i]::complex(numeric) then
                     print("C",[theta_init*Efun[i], z = 0 .. d]);
                  end if;
               end if;
               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 solve A*x=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):

##
             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:=CodeTools:-Usage([doCalc(-0.06, 0.001)]):
evalf(ans[9]);  ## H value
evalf(ans[10]);  ## H bar

memory used=451.31MiB, alloc change=336.20MiB, cpu time=2.84s, real time=3.06s, gc time=250.00ms

with(plots):

d:= 0.2e-3:

testproc := proc(j, u, k) option remember;
   local calcvals;
   if not [j,u,k]::[numeric,numeric,posint] then
      return 'procname'(args);
   end if;
   calcvals:=doCalc(j,u);
   evalf(calcvals[k]);
end proc:

NN := 8;

(gridji,rngj,rngi):=[51,101],0..3,-2.0..-0.0;
PP[gridji,4]:=plot3d(max(-1e2,testproc(i,j,4)), j=rngj, i=rngi, grid=gridji):
PP[gridji,7]:= plot3d(Im(testproc(i,j,7)), j=rngj, i=rngi, grid=gridji):

## Surface plot using 3d
PP[gridji,4];

## This is the pmin values
data[gridji,4]:=op([1,3],PP[gridji,4]):

## Convert Array to Matrix to allow the use listcontplot
ConMatrix := Matrix(data[gridji,4]);

## Assign the min and max value of pmin to minv and maxv
#(minv,maxv):=[min,max](data[gridji,4])[];

## GRB color and the contours
(color1,color2):="Red","Yellow":
conts:=[minv, seq(minv+(i-1)*(maxv-minv)/(NN+2-1),i=2..NN+1), maxv];

colorlist:=ColorTools:-Color~(((c1,c2,N)->[seq(c1+(i-1)*(c2-c1)/(N-1),
                                               i=1..N)])([ColorTools:-Color(color2)[]],
                                                         [ColorTools:-Color(color1)[]],
                                                         nops(conts))):

##

## Use the Conmatrix and the conts to plot contours and transform the axes
subsindets(plots:-listcontplot(ConMatrix, contours=conts,filledregions,thickness=0),
           specfunc(anything,THICKNESS),u->THICKNESS(0.2)):
LCP[gridji,4]:=plottools:-transform((x,y)->[lhs(rngj)+(rhs(rngj)-lhs(rngj))*(x-1)/(gridji[1]-1),
                                            lhs(rngi)+(rhs(rngi)-lhs(rngi))*(y-1)/(gridji[2]-1)])(%):
display(LCP[gridji,4],
        seq(plot(-2.0,1..1,legend=evalf[4](conts[i]),thickness=15,
                 color=colorlist[nops(conts)-i+1]),
            i=1..nops(conts)),
        size=[500,400],legendstyle=[location=right], axes=boxed, labels = ["u", "xi"], labelfont = ["TimesNewRoman", 14], axesfont = [14, 14]);

## This deals with the Imaginary part

## Surface plot using plot3d for Imginary part
PP[gridji,7]:

##

data[gridji,7]:=op([1,3],PP[gridji,7]):

## Convert Array to matrix
AMatrix := Matrix(data[gridji,7]);

(minv,maxv):=[min,max](data[gridji,7])[];

(color1,color2):="Red","Yellow":
conts1:=[minv, seq(minv+(i-1)*(maxv-minv)/(NN+2-1),i=2..NN+1), maxv]/~1.5;

colorlist:=ColorTools:-Color~(((c1,c2,N)->[seq(c1+(i-1)*(c2-c1)/(N-1),
                                               i=1..N)])([ColorTools:-Color(color2)[]],
                                                         [ColorTools:-Color(color1)[]],
                                                         nops(conts1))):

subsindets(plots:-listcontplot(AMatrix,
                               contours=conts1-~1e-9,filledregions,thickness=0),
           specfunc(anything,THICKNESS),u->THICKNESS(0.2)):
LCP[gridji,7]:=plottools:-transform((x,y)->[lhs(rngj)+(rhs(rngj)-lhs(rngj))*(x-1)/(gridji[1]-1),
                                            lhs(rngi)+(rhs(rngi)-lhs(rngi))*(y-1)/(gridji[2]-1)])(%):
display(LCP[gridji,7],
        seq(plot(-2.0,1..1,legend=evalf[4](conts1[i]),thickness=15,
                 color=colorlist[nops(conts)-i+1]),
            i=1..nops(conts1)),
        size=[500,400],legendstyle=[location=right], axes=boxed, labels = ["u", "xi"], labelfont = ["TimesNewRoman", 14], axesfont = [14, 14]);