restart;

  A := Matrix([[2*t, 2, 3], [4, 5*t, 6], [7, 8*t, 9]]);
  B := Matrix([t, 2*t, 3*t])^%T;
  test1:=proc(n)
              local tau, s,M,V;
              global t;
              s:=1:
              for t from 1 to n do
                  M:=A+A^(-1);
                  V[s]:=(M.B):
                  s:=s+1:
              end do;
              V;
         end proc:

test2:=proc(n)
              local s,M,V;
              global t;
              s:=1:
              for t from 1 to n do
                  M:=A+A^(-1);
                  V[s]:=(M.B):
                  s:=s+1:
              end do;
              eval(V);
         end proc:

  R:=test1(4);
  R[1];

  R:=test2(4);

restart:

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:= 30:

                        for i from 1 to m do
                        w := 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];

  restart:
  with(plots): with(geometry):
  _EnvHorizontalName := 'x':
  _EnvVerticalName := 'y':
  a := 7:
  point(E, 0, a*sqrt(3)/2):
  point(B, -a/2, 0):
  point(C, a/2, 0):
  point(o, 0, a*sqrt(3)/6):
  point(A, 0, a/2):
  point(H, 0, 0):
  R := (3-sqrt(3))*sqrt(2)*a/12:
  point(J, 0, a*sqrt(3)/6 - R):
  line(L1, -(7*x*sqrt(3))/6 + (7*y)/2 - (49*sqrt(3))/12 = 0):
  reflection(J1, J, L1):
  coordJ1:=coordinates(J1):
  
  line(L2, -(7*x*sqrt(3))/6 - (7*y)/2 + (49*sqrt(3))/12 = 0):
  reflection(J2, J, L2):coordJ2:=coordinates(J2):
  triangle(Tr1, [E, B, C]):
  triangle(Tr2, [A, B, C]):
  StretchRotation(E1, E, B, Pi/4, clockwise, sqrt(2)/2):
  coordinates(E1):
  StretchRotation(E2, E, C, Pi/4, counterclockwise, sqrt(2)/2):
  coordinates(E2):
  triangle(Tr3, [E, B, E1]):
  triangle(Tr4, [E, C, E2]):
  triangle(Tr5, [B, C, J]):
  circle(cir, [point(P1,[0,a*sqrt(3)/6]), R]):
  poly1 := Matrix( [ [0, a*sqrt(3)/2],
                     coordJ1,
                    [0, a/2],
                    coordJ2
                  ],
                  datatype = float
                ):
  pol1 := polygonplot(poly1, colour = green, transparency = 0.7, gridlines):
  pol2:= plottools:-rotate(pol1, 5*Pi/3, coordinates(o)):
  display( [pol1, pol2] );
(*  poly2 := Matrix([ coordinates(J),coordinates(E1), coordJ2,|-a/2,0]],
                  datatype = float
                ):
    pol2 := polygonplot(poly2, colour = green, transparency = 0.7, gridlines):
*)
  tex := textplot([[0, a*sqrt(3)/2, "E"],[0,a*sqrt(3)/6 - R,"J"],
                    [-7/4 + (7*sqrt(3))/4, -7/4 + (7*sqrt(3))/4,"E1"],
                    [7/4 - (7*sqrt(3))/4, -7/4 + (7*sqrt(3))/4,"E2"],
                    [-7/4 + (((7*sqrt(3))/6 - 7*(3 - sqrt(3))*sqrt(2)/12)*sqrt(3))/2,
                    (7*sqrt(3))/6 + 7*(3 - sqrt(3))*   sqrt(2)/24,"J1"],
                    [7/4 - (((7*sqrt(3))/6 - 7*(3 - sqrt(3))*sqrt(2)/12)*sqrt(3))/2,
                    (7*sqrt(3))/6 + 7*(3 - sqrt(3))*   sqrt(2)/24,"J2"]], 'align' = {'above', 'right'}):
  display( [ draw( [ Tr1(color = cyan),
                     Tr3(color = green),
                     cir(color=blue),
                     Tr2(color = red),
                     Tr4(color = grey),
                     Tr1(color=blue)],'view' = [-5 .. 5, 0 .. 7
                   ],
                   axes = normal,
                   scaling = constrained,
                   size=[800,800]
                   ),
            tex
           ]
          );
  line(L1, [B, o]):
  Equation(L1):
  line(L1, -(7*x*sqrt(3))/6 + (7*y)/2 - (49*sqrt(3))/12 = 0):
  reflection(J1, J, L1):coordinates(J1):
  triangle(Tr6, [B, J1, E]):
  line(L2, [C, o]):
  Equation(L2):
  line(L2, -(7*x*sqrt(3))/6 - (7*y)/2 + (49*sqrt(3))/12 = 0):
  reflection(J2, J, L2):coordinates(J2):
  triangle(Tr7, [C, J2, E]):
  triangle(T1, [E, J1, A]):
  triangle(T2, [E, C, E2]):
  triangle(T3, [B, H, J]):
  triangle(T4, [C, H, J]):

  display( [ draw( [ cir(color = orange, filled = true, transparency = 0.1),
                     Tr6(color = blue, filled = true, transparency = 0.2),
                     Tr5(color = blue, filled = true, transparency = 0.2),
                     Tr7(color = blue, filled = true, transparency = 0.2)
                     
                   ],
                   axes = none,
                   scaling = constrained
                 ),
             pol1
           ]
         );

  restart:
  with(plots): with(geometry):
  _EnvHorizontalName := 'x':
  _EnvVerticalName := 'y':
  a := 7:
  point(E, 0, a*sqrt(3)/2):
  point(B, -a/2, 0):
  point(C, a/2, 0):
  point(o, 0, a*sqrt(3)/6):
  point(A, 0, a/2):
  point(H, 0, 0):
  R := (3-sqrt(3))*sqrt(2)*a/12:
  point(J, 0, a*sqrt(3)/6 - R):
  line(L1, -(7*x*sqrt(3))/6 + (7*y)/2 - (49*sqrt(3))/12 = 0):
  reflection(J1, J, L1):
  coordJ1:=coordinates(J1):
  
  line(L2, -(7*x*sqrt(3))/6 - (7*y)/2 + (49*sqrt(3))/12 = 0):
  reflection(J2, J, L2):coordJ2:=coordinates(J2):
  triangle(Tr1, [E, B, C]):
  triangle(Tr2, [A, B, C]):
  StretchRotation(E1, E, B, Pi/4, clockwise, sqrt(2)/2):
  coordinates(E1):
  StretchRotation(E2, E, C, Pi/4, counterclockwise, sqrt(2)/2):
  coordinates(E2):
  triangle(Tr3, [E, B, E1]):
  triangle(Tr4, [E, C, E2]):
  triangle(Tr5, [B, C, J]):
  circle(cir, [point(P1,[0,a*sqrt(3)/6]), R]):
  poly1 := Matrix( [ [0, a*sqrt(3)/2],
                     coordJ1,
                    [0, a/2],
                    coordJ2
                  ],
                  datatype = float
                ):
  pol1 := polygonplot(poly1, colour = green, transparency = 0.7, gridlines):
(*  poly2 := Matrix([ coordinates(J),coordinates(E1), coordJ2,|-a/2,0]],
                  datatype = float
                ):
    pol2 := polygonplot(poly2, colour = green, transparency = 0.7, gridlines):
*)
  tex := textplot([[0, a*sqrt(3)/2, "E"],[0,a*sqrt(3)/6 - R,"J"],
                    [-7/4 + (7*sqrt(3))/4, -7/4 + (7*sqrt(3))/4,"E1"],
                    [7/4 - (7*sqrt(3))/4, -7/4 + (7*sqrt(3))/4,"E2"],
                    [-7/4 + (((7*sqrt(3))/6 - 7*(3 - sqrt(3))*sqrt(2)/12)*sqrt(3))/2,
                    (7*sqrt(3))/6 + 7*(3 - sqrt(3))*   sqrt(2)/24,"J1"],
                    [7/4 - (((7*sqrt(3))/6 - 7*(3 - sqrt(3))*sqrt(2)/12)*sqrt(3))/2,
                    (7*sqrt(3))/6 + 7*(3 - sqrt(3))*   sqrt(2)/24,"J2"]], 'align' = {'above', 'right'}):
  display( [ draw( [ Tr1(color = cyan),
                     Tr3(color = green),
                     cir(color=blue),
                     Tr2(color = red),
                     Tr4(color = grey),
                     Tr1(color=blue)],'view' = [-5 .. 5, 0 .. 7
                   ],
                   axes = normal,
                   scaling = constrained,
                   size=[800,800]
                   ),
            tex
           ]
          );
  line(L1, [B, o]):
  Equation(L1):
  line(L1, -(7*x*sqrt(3))/6 + (7*y)/2 - (49*sqrt(3))/12 = 0):
  reflection(J1, J, L1):coordinates(J1):
  triangle(Tr6, [B, J1, E]):
  line(L2, [C, o]):
  Equation(L2):
  line(L2, -(7*x*sqrt(3))/6 - (7*y)/2 + (49*sqrt(3))/12 = 0):
  reflection(J2, J, L2):coordinates(J2):
  triangle(Tr7, [C, J2, E]):
  triangle(T1, [E, J1, A]):
  triangle(T2, [E, C, E2]):
  triangle(T3, [B, H, J]):
  triangle(T4, [C, H, J]):

  display( [ draw( [ cir(color = orange, filled = true, transparency = 0.1),
                     Tr6(color = blue, filled = true, transparency = 0.2),
                     Tr5(color = blue, filled = true, transparency = 0.2),
                     Tr7(color = blue, filled = true, transparency = 0.2)
                     
                   ],
                   axes = none,
                   scaling = constrained
                 ),
             pol1
           ]
         );

  restart:
  N := 2.7:
#
# Data as lists
#
  Hb := [0.076, 0.083]:
  k := [0.00003400566801, 0.00003424620533]:
  P50a := [20.78938475, 21.39546041]:
  P50v := [21.20711722, 22.06793197]:
  nu := [0.02042461957, 0.02120393111]:
  Df := [0.00001617321837, 0.00001607066092]:
  P__baro := [759.062, 759.062]:
  PaCor := [94.82734101, 90.40928915]:
  PvCor := [35.32630403, 35.55779803]:
 

#
# My original calculation, using seq() commands to
# step through lists
#
  f := [ seq( 1.34*N*Hb[j]*(p/((P50a[j] + P50v[j])/2))^(N - 1)/((P50a[j] + P50v[j])/2*k[j]*(1 + (p/((P50a[j] + P50v[j])/2))^N)^2),
              j = 1 .. 2
            )
       ]:
  result := [seq( Por*phi*(nu[j]/Df[j])^(2/3)*int( (1 + f[j])^(2/3)/(P__baro[j] - p),
                                                    p = PvCor[j] .. PaCor[j]
                                                 )/(4*(1 - Por)*BP__length),
                  j = 1 .. 2
                )
            ];

#
# Carl Love's code - define a procedure for the calculation and
# apply it elementwise (works for lists as well as vectors)
#
  MyProc:= proc(N, Hb, k, P50a, P50v, nu, Df, P__baro, PaCor, PvCor)
                local p, P:= (2*p/(P50a+P50v))^N;
                Por*phi/4/(1-Por)/BP__length*(nu/Df)^(2/3) * int( (1+1.34*Hb*N*P/k/p/(1+P)^2)^(2/3)/(P__baro - p),
                                                                   p= PvCor..PaCor
                                                                )
           end proc:
  MyProc~(N, Hb, k, P50a, P50v, nu, Df, P__baro, PaCor, PvCor);

restart:

Procedure

Run the calculation for supplied value of 'xi'

r(q) graph

  restart:
  with(plots):
  with(plottools):
  N:=10:
  t_upper:=evalf(Pi):
  s_step:=0.01:
   
  PDE:=diff(u(y, t), t) = diff(u(y, t), y, y);
  ICBC:= [u(y,0) = 0, u(0,t) = cos(t), u(N, t) = 0];
#
# Note that boundary conditions are contradictory. The
# first requires u(0,0)=0 and the second requires
# u(0,0)=cos(0)=1. Not sure what Maple is going to do
# with this

#
# Set up plots of the boundary conditions so that
# these can be included in a general plot of the
# solution surface - just as a check!!!
#
  f1:= plottools:-transform((x, y) -> [0, x, y]):
  f2:= plottools:-transform((x, y) -> [x, 0, y]):
  q1:=plot( cos(x), x=0..t_upper):
  q2:=plot( 0,      x=0..t_upper):
  q3:=plot( 0,      x=0..N):
  pbc:=display( [ f1(q1), f1(q2), f2(q3) ],
                color=green,
                thickness=4
              );

#
# Solve system using Maple defaults for spacestep and timestep
# and plot with boundary conditions. Note that spacestep will
# default to (N-0)/20, ie 0.5 and timestep will default to this
# value of spacestep
#
# Looks like Maple is *trying* to use the boundary condition
# u(0, t)=cos(t) along y=0 - except close to t=0, where it is
# using u(0,0)=0
#
  sol1:= pdsolve(PDE,ICBC,numeric):
  display( [ sol1:-plot3d( u(y, t),
                           y=0..N,
                           t=0..t_upper,
                           color=red,
                           style=surface
                         ),
              pbc
            ],
            title="Maple defaults, check region along y=0, t=0..1",
            titlefont=[times, roman, 20]
          );

#
# So check the returned value for t=0, y=0
#
# 1. from the boundary condition  u(y, 0)=0 the
#    answer for u(0,0) should return 0
# 2. from the boundary condition  u(y, 0)=cos(t)
#    the answer should be 1, when t=0
#
  uVal := sol1:-value(u(y,t), t=0):
  eval(u(y,t), uVal(0));
#
# So we seem to be obeying the u(y, 0)=0 condition
#
  

#
# In the above the boundary condition u(0,t) = cos(t), can get a
# better fit along y=0, by reducing the spacestep (which will
# reduce the timestep, and the reduced timestep will improve
# the fit along u(0,t) = cos(t)
#
  sol2:= pdsolve(PDE,ICBC,numeric, spacestep=s_step):
  display( [ sol2:-plot3d( u(y, t),
                           y=0..N,
                           t=0..t_upper,
                           color=blue,
                           style=surface
                         ),
              pbc
            ],
            title="spacestep=timestep=0.1, compare region along y=0, t=0..1 with above",
            titlefont=[times, roman, 20]
          );
 

#
# So check the returned value for t=0, y=0
#
# 1. from the boundary condition  u(y, 0)=0 the
#    answer for u(0,0) should return 0
# 2. from the boundary condition  u(y, 0)=cos(t)
#    the answer should be 1, when t=0
#
  uVal := sol2:-value(u(y,t), t=0):
  eval(u(y,t), uVal(0));
#
# So we seem to be obeying the u(y, 0)=0 condition
#

#
# Plot the both solutions along the line y=0, for comparison with
# the boundary condition u(0,t)=cos(t), Both solutions are very
# wrong at t=0, although sol2 (with reduced spacestep and timestep)
# is much better thereafter. Note that by default sol1, will have
# a (default) timestep of 0.5, whereas for sol2 it is set by s_step
# (defined as 0.01 above)
#
  display( [ sol1:-plot(y=0, t=0..t_upper, color=red),
             sol2:-plot(y=0, t=0..t_upper, color=blue),
             plot( cos(t), t=0..t_upper, color=green)
           ]
         );

restart;

N:=10;

PDE:=diff(u(y, t), t) = diff(u(y, t), y, y);

ICBC:= {u(y,0) = 0, u(0,t) = cos(t), u(N, t) = 0};

dy := N/100;  # space step
dt := dy^2;   # time step

sol1:=pdsolve(PDE,ICBC,numeric,spacestep=dy, timestep=dt):

sol1:-plot3d(t=0..Pi);

vals := sol1:-value(u(y,t), output=listprocedure);

U := eval(u(y,t), vals);

U(1,0.5);

#
# Solution at (0,0) - please explain how this statisfies
# the BC u(0,t) = cos(t)
#
  U(0,0);
#
# Demonstrate the discrepancy along the line x=0
#
  plots:-display( [ plot( U(0,t), t=0..100*dt, color=red),
                    plot( cos(t), t=0..100*dt, color=blue)
                  ]
                );
#
# Check along the line t=0. Note that the boundary condition
# u(0,t) = cos(t), requires U(0,0)=1
#
  plot( U(x,0),x=0..N, color=red, axes=boxed)

-D[1](U)(1,0.5);

plot3d(-D[1](U)(y,t), y=0..N, t=0..Pi);