Please, how do I use surface or contour to plot discontinous function?

z = f(x,y) := x-.8193318913*sin(x)*cos(x)/(cos(x)^2+sinh(y)^2)-(0.2931044702e-2*(0.7500000000e-3*x^2-0.7500000000e-3*y^2+0.1622336517e-1))*sin(x)*cos(x)/(((0.7500000000e-3*x^2-0.7500000000e-3*y^2+0.1622336517e-1)^2+0.2250000000e-5*x^2*y^2)*(cos(x)^2+sinh(y)^2))-0.4396567053e-5*x*y*sinh(y)*cosh(y)/(((0.7500000000e-3*x^2-0.7500000000e-3*y^2+0.1622336517e-1)^2+0.2250000000e-5*x^2*y^2)*(cos(x)^2+sinh(y)^2)):

Hi, 

Please how I do plot: theta := -(65.7014900075861*(cos(-4.536529763+45365.29764*z)+.1749541674))*exp(-1.603200636*t) for z=0..d.

I tried this: 

display(plot([seq(subs(t = i, theta), i = [seq(0.1*ii, ii = 1 .. 7)])], z = 0 .. d));

But I want this sequence: [seq(subs(t = i, theta), i = [seq(0.1*ii, ii = 1 .. 7)])] to include Pi/2 such that t_0 = Pi/2 in the plot. I.e., I want the Pi/2 to be the initial value in the sequence.

The matrix data[gridji,4] in my code has float undefined. Please, how do I handle it?

I want to  return {x = 0., y = 1.158748796 + 0. I}  as  {x = 0., y = 1.158748796 }.  The solution is coming from:

soln3:= fsolve({b1, b2}, {x = 0 .. infinity, y = 0 .. infinity}); 

and the second solution is coming from:
 soln4:= fsolve({b1, b2}, {x = -infinity ..0, y = -infinity .. 0}); 

See my code below

restart:

Procedure

doCalc:= proc( xi )

                 # 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.0001,
                       eta__1:= 0.240e-1,
                       alpha:= 1-alpha__3^2/(gamma__1*eta__1), 
                       c:= alpha__3*xi*alpha/(eta__1*(4*k__1*q^2/d^2-alpha__3*xi/eta__1)),
                       theta_init:= theta0*sin(Pi*z/d),
                       n:= 30,
                       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, soln4, Imagroot1, Imagroot2;

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

# Compute the rest of the Roots
  soln3:= fsolve({b1, b2}, {x = 0 .. infinity, y = 0 .. 10}); 
  soln4:= fsolve({b1, b2}, {x = -infinity ..0, y = -infinity .. 0}); 
  Imagroot1:=subs(soln3, I*y); 
  Imagroot2:= subs(soln4, I*y); 
  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):

## Return all the plots
            return qq, Imagroot1,Imagroot2, p, soln3, soln4;
  end proc:

ans:=[doCalc(0.06)]:
ans[5];
                {x = 0., y = 1.158748796 + 0. I}
ans[6];
                   {x = 0., y = -1.158748796}
 

I'm using a procedure but the out is not returned. Please why?

restart:

Procedure

doCalc:= proc( xi )

                 # 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.0001,
                       eta__1:= 0.240e-1,
                       alpha:= 1-alpha__3^2/(gamma__1*eta__1), 
                       c:= alpha__3*xi*alpha/(eta__1*(4*k__1*q^2/d^2-alpha__3*xi/eta__1)),
                       theta_init:= theta0*sin(Pi*z/d),
                       n:= 30,
                       g, f, b1, b2, qstar, OddAsymptotes, ModifiedOddAsym,
                       qstarTemporary, indexOfqstar2, qstar2, AreThereComplexRoots,
                       soln1, soln2, qcomplex1, qcomplex2, gg, qq, m, pp, j, i,
                       AllAsymptotes, w, Efun, b, aa, F, A, B, Ainv, r, theta_sol, v, Vfun, v_sol,minp,nstar,
                       soln3, soln4, Imagroot1, Imagroot2;

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

                       qq := [2.106333379+.6286420119*I, 2.106333379-.6286420119*I, 4.654463885, 7.843624703, 10.99193295,14.13546782, 17.27782732, 20.41978346, 23.56157073, 26.70327712, 29.84494078, 32.98658013,         36.12820481, 39.26982019, 42.41142944, 45.55303453, 48.69463669, 51.83623675, 54.97783528,                     58.11943264, 61.26102914, 64.40262495, 67.54422024, 70.68581510, 73.82740963,                                 76.96900389, 80.11059792, 83.25219177, 86.39378546, 89.53537903, 92.67697249];
                        m:= numelems(qq):

                        for i from 1 to m do 
                        w[i] := gamma__1*alpha/(4*k__1*qq[i]^2/d^2-alpha__3*xi/eta__1);
                        end do;

## Return all the plots
  return w;
  end proc:

ans:=[doCalc(0.06)]:
ans[1];