DrugTrial:= <
   20, 11, 19;
   4,  4,  17
>:

infolevel[ContingencyTableTests]:= 1:

ContingencyTableTests(DrugTrial, method= Pearson):  % = evalf(%);

#Compare with:
Statistics:-ChiSquareIndependenceTest(DrugTrial);

infolevel[ContingencyTableTests]:= 0:
ContingencyTableTests(DrugTrial, method= Yates):  % = evalf(%);

ContingencyTableTests(DrugTrial, method= G):  % = evalf(%);

CodeTools:-Usage(ContingencyTableTests(DrugTrial, method= Fisher)):  % = evalf(%);

can greatly facilitate the organization of one's work,

can encapsulate the setting of parameter values,

can allow one to work with symbolic parameters.

restart:

KandR:= module()
local
   a, b, c, e, #parameters

   #procedure that lets user set parameter values:
   ModuleApply:= proc({
       a::algebraic:= KandR:-a, b::algebraic:= KandR:-b,
       c::algebraic:= KandR:-c, e::algebraic:= KandR:-e
   })
   local k;
      for k to _noptions do thismodule[lhs(_options[k])]:= rhs(_options[k]) od;
      return
   end proc,

   Extinct:= (x::realcons)-> `if`(x < 0.5, 0, x) #force small, insignificant values to 0
;
export
   #Procedure that does one symbolic iteration
   #(Note that this procedure uses Vector input and output.)
   iter_symb:= KR-> KR *~ <-b*KR[2]+a+1, b*e*KR[1]-c+1>, 

   #Such simple treatment as above is only possible for autonomous
   #recurrences.

  
   iter_num:= Extinct~@iter_symb #one numeric iteration
;
end module:

#The following expression is the discrete equivalent of the derivative (or gradient).
#It represents the change over one time period.
P:= <K,R>:  
OneStep:= KandR:-iter_symb(P) - P

#An equilibrium occurs when the gradient is 0.
Eq:= <K__e, R__e>:
Eqs:= solve({seq(eval(OneStep=~ 0, [seq(P=~ Eq)]))}, [seq(Eq)]);

#We're only interested here in nonzero solutions.
EqSol:= remove(S-> 0 in rhs~(S), Eqs)[];

#Set parameters:
KandR(a= 0.15, b= 0.01, c= 0.02, e= 0.01);

#Show idempotency at equilibrium:
use KandR in Eq0:= eval(Eq, EqSol); print(Eq0 = iter_num(Eq0)) end use:

#procedure that fills a Matrix with computed values of a 1st-order recurrence.
#(A more-efficient method than this can be used for linear recurrences.)
#This procedure has no dependence on the module.
Iterate:= proc(n::nonnegint, iter, init::Vector[column])
local M:= Matrix((n+1, numelems(init)), init^+, datatype= hfloat), i;
   for i to n do M[i+1,..]:= iter(M[i,..]) od;
   M
end proc:

dev:= <7,2>:
use KandR in KR:= Iterate(1000, iter_num, Eq0 + dev) end use:

plot(
   [
       KR, #trajectory of population
       KR[[1,1],..], #1st point
       KR[-[1,1],..], #last point,
       <Eq0|Eq0>^+, #equilibrium
       #every 100th point (helps show time scale):
       KR[100*[$1..iquo(numelems(KR[..,1]), 100)-1], ..]
   ],
   #This group of options are all lists, each element of which corresponds
   #to one of the above components of the plot:
   style= [line, point$4],
   symbol= [solidcircle$4, soliddiamond],
   symbolsize= [18$4, 12],
   color= [black, green, red, brown, blue],
   thickness= [0$5],
   legend= [`pop.`, init, final, equilibrium, `100 periods`],

   #This group of options are lists, each element of which corresponds to one
   #coordinate axis (horizontal, then vertical).
   view= [0..max(KR[..,1]), 0..max(KR(..,2))],
   labels= [rabbits, foxes],
   labeldirections= [horizontal,vertical],
   size= [700,700], #measured in pixels

   #options applied to whole plot:
   labelfont= [TIMES, BOLDITALIC, 14],
   title= "Population of foxes and rabbits over time" "\n", titlefont= [TIMES,16],
   caption=
      "\n" "Choosing an initial point near the equilibrium causes"
      "\n" "outward spiraling divergence." "\n",
   gridlines= false
);
 

plots:-fieldplot(
   rtable_eval(OneStep),
   K= 0..max(KR[..,1]),  R= 0..max(KR[..,2]), grid= [16,16],

   #arrow-specific options:
   anchor= tail, fieldstrength= log, arrows= slim, color= "DarkGreen",

   #other options (same as any 2D plot):
   labels= [rabbits, foxes], labeldirections= [horizontal,vertical],
   labelfont= [TIMES, BOLDITALIC, 14],
   title= "One-step population changes from any point" "\n", titlefont= [TIMES,16],
   caption= "\n" "All trajectories spiral outward from the equilibrium." "\n",
   size= [700,700],
   gridlines= false
);

restart:

Digits:= trunc(evalhf(Digits)); #generally a very efficient setting

#ordinary differential equations:
ODEs:= [
   #Eq 8:
   diff(f(eta),eta$4)
   - S*(diff(f(eta),eta$3)*(eta-f(eta)) + diff(f(eta),eta$2)*(3+diff(f(eta),eta)))
   - Ha^2*diff(f(eta),eta$2),

   #Eq 9:
   (1+4/3*Rd)*diff(theta(eta),eta$2)
   + Pr*(S*diff(theta(eta),eta)*(f(eta)-eta) + Ec*diff(f(eta),eta$2)^2)
   + Nb*diff(phi(eta),eta)*diff(theta(eta),eta) + Nt*diff(theta(eta),eta)^2,

   #Eq 10:
   diff(phi(eta),eta$2) + S*Sc*diff(phi(eta),eta)*(f(eta)-eta)
   + Nt/Nb*diff(theta(eta),eta$2)

   #All these ODEs are implicitly equated to 0.
]:

<ODEs[]>; #Display the ODEs.
   

#lower and upper boundary values of the independent variable:
(LB,UB):= (0,1): #Eq 11

#boundary conditions:
BCs:= [
   #Eq 11:
   ([f, (D@@2)(f), D(theta), D(phi)](LB) =~ [0, 0, 0, 0])[],
   ([f, D(f), theta, phi](UB) =~ [1, 0, 1, 1])[]
];

#parameters (input values to BVP system) and their default values:
Params:= Record(
   #Eq 12:
   S=  0.5,  #Squeeze number
   Pr= 10.,  #Prandtl number
   Ec= 0.1,  #Eckert number
   Sc= 0.5,  #Schmidt number
   Ha= 2.,   #Hartman number
   Nb= 0.15, #Brownian motion parameter
   Nt= 0.15, #thermophoretic parameter
   Rd= 4.5   #radiation parameter
):
   

#Named special computed boundary values (output values of the system): Note that each adds
#exactly one parameter and one boundary condition to the system.

NBVs:= [   
   -(D@@2)(f)(1) = `C*__f`,  #Eq 13: Skin friction coefficient
   -D(theta)(1) = `Nu*`      #Eq 14: Nusselt number
]:

#just rewriting the above names in a more-convenient form:
Nu:= `Nu*`:
Cf:= `C*__f`:

#Appliable module that calls dsolve for each
#particular configuration of parameter values. (Maple's numeric BVPs require a separate
#dsolve invocation for each specific set of parameter values (unlike IVPs, which can be
#parameterized into a single dsolve invocation).)
#
#For every dsolve call through this module, the NBVs are evaluated and saved.

Solve:= module()
local
   nbvs_rhs:= rhs~(:-NBVs), #just the names
   Sol, #numeric dsolve BVP solution of any 'output' form
   ModuleApply:= subs(
      _Sys= {:-ODEs[], :-BCs[], :-NBVs[]},
      proc({
         #These are just default parameter values: They can be changed when the
         #procedure is called.
         S::realcons:=  Params:-S,
         Pr::realcons:= Params:-Pr,
         Ec::realcons:= Params:-Ec,
         Sc::realcons:= Params:-Sc,
         Ha::realcons:= Params:-Ha,
         Nb::realcons:= Params:-Nb,
         Nt::realcons:= Params:-Nt,
         Rd::realcons:= Params:-Rd
         #By extracting the default values from record Params, it is possible to change
         #them by changing the record.
      })
         Sol:= dsolve(_Sys, _rest, numeric);
         AccumData(Sol, {_options});
         Sol
      end proc
   ),
   AccumData:= proc(
      Sol::{Matrix, procedure, list({name, function}= procedure)},
      params::set(name= realcons)
   )
   local n, nbvs;
      if Sol::Matrix then
         nbvs:= seq(n = Sol[2,1][1,Pos(n)], n= nbvs_rhs)
      else
         nbvs:= (nbvs_rhs =~ eval(nbvs_rhs, Sol(:-LB)))[]
      fi;
      SavedData[params]:= Record[packed](params[], nbvs)
   end proc,
   ModuleLoad:= eval(Init)
;
export
   SavedData, #table of Records
   Pos, #Matrix column indices of nbvs
   Init:= proc()
      Pos:= proc(n::name) option remember; local p; member(n, Sol[1,1], 'p'); p end proc;
      SavedData:= table();
      return
   end proc   
;
   ModuleLoad()
end module:
 

colseq:= [red, green, blue, brown]:

#parameter values that remain fixed for the entire set of plots:
Pc:= Sc= 0.5, Nb= 0.1, Ec= 0.1, Pr= 10:
 

#parameter values that remain fixed with each of the four plots::
Ps:= [
   [S= 0.5, Ha= 2, Nt= 0.1],
   [Ha= 2, Rd= 0.5, Nt= 0.1],
   [S= 0.5, Rd= 0.5, Nt= 0.1],
   [S= 0.5, Ha= 2, Rd= 0.5]
]:

#parameter value for each curve
Pv:= [
   Rd= [0.5, 2, 4, 8],
   S=  [2, 4, 8, 12],
   Ha= [2, 4, 8, 12],
   Nt= [0.01, 0.05, 0.1, 0.2]
]:
      

for i to nops(Ps) do
   plots:-display(
      [seq(
         plots:-odeplot(
            Solve(lhs(Pv[i])= rhs(Pv[i])[j], Ps[i][], Pc),
            [eta, theta(eta)], 'color'= colseq[j], 'legend'= [lhs(Pv[i])= rhs(Pv[i])[j]]
         ),
         j= 1..nops(rhs(Pv[i]))
      )],
      'axes'= 'boxed', 'gridlines'= false,
      'labelfont'= ['TIMES', 'BOLDOBLIQUE', 16],
      'caption'= nprintf(
         cat("\n%a = %4.2f, "$nops(Ps[i])-1, "%a = %4.2f\n\n"), (lhs,rhs)~(Ps[i])[]
      ),
      'captionfont'= ['TIMES', 16]
   )
od;

#procedure that generates a 2-D plot of an expression versus a varied parameter:
ParamPlot2d:= proc(
   Y::{procedure, `module`}, #procedure that extracts y-value from Solve's dsolve solution
   X::name= range(realcons), #x-axis-parameter range
   FP::list(name= realcons), #fixed values of other parameters
   {
      eta::realcons:= :-LB, #independent variable value
      dsolveopts::list({name, name= anything}):= []  
   }
)
   plot(
      x-> Y(
         Solve(
            lhs(X)= x, FP[],
            #Default dsolve options can be changed by setting 'dsolveopts':
            'abserr'= 0.5e-7, 'interpolant'= false, 'output'= Array([eta]), dsolveopts[]
         )
      ),
      rhs(X),
      #Default plot options can be changed by putting the option in the ParamPlot2d call:
      'numpoints'= 25, 'axes'= 'boxed', 'gridlines'= false,
      'labelfont'= ['TIMES', 'BOLDOBLIQUE', 16],
      'caption'= nprintf(cat("%a = %4.2f, "$nops(FP)-1, "%a = %4.2f"), (lhs,rhs)~(FP)[]),
      'captionfont'= ['TIMES', 16],
      _rest
   )
end proc:

#procedure that extracts Nusselt number from dsolve solution:
GetNu:= (Sol::Matrix)-> Sol[2,1][1, Solve:-Pos(:-Nu)]:

RD:= [0.5, 2, 4]:
plots:-display(
   [seq(
      ParamPlot2d(
         GetNu, Ha= 2..8, [S= 0.5],
         dsolveopts= [Rd= RD[k], Sc= 0.5, Nb= 0.1, Ec= 0.1, Nt= 0.1, Pr= 10],
         'legend'= [Rd= RD[k]], 'color'= colseq[k], 'labels'= [Ha, Nu]
      ),
      k= 1..nops(RD)
   )]
);

SV:= [0.5, 2, 4]:
plots:-display(
   [seq(
      ParamPlot2d(
         GetNu, Ha= 2..8, [Rd= 0.5],
         dsolveopts= [S= SV[k], Sc= 0.5, Nb= 0.1, Ec= 0.1, Nt= 0.1, Pr= 10],
         'legend'= [S= SV[k]], 'color'= colseq[k], 'labels'= [Ha, Nu]
      ),
      k= 1..nops(SV)
   )]
)

#procedure that generates 3-D plots (dropped-shadow contour + surface) of an expression
#versus two varied parameters:

ParamPlot3d:= proc(
   Z::{procedure, `module`}, #procedure that extracts z-value from Solve's dsolve solution
   X::name= range(realcons), #x-axis-parameter range
   Y::name= range(realcons), #y-axis-parameter range
   FP::list(name= realcons), #fixed values of other parameters
   {
      #fraction of empty space above and below plot (larger "below"
      #value improves view of dropped-shadow contourplot):
      zmargin::[realcons,realcons]:= [.05,0.15],
      eta::realcons:= :-LB, #independent variable value
      dsolveopts::list({name, name= anything}):= [],
      contouropts::list({name, name= anything}):= [],
      surfaceopts::list({name, name= anything}):=[]    
   }
)
local
   LX:= lhs(X), RX:= rhs(X), LY:= lhs(Y), RY:= rhs(Y),
   Zremember:= proc(x,y)
   option remember; #Used because 'grid' should be the same for both plots.
      Z(
         Solve(
            LX= x, LY= y, FP[],
            #Default dsolve options can be changed by setting 'dsolveopts':
            'abserr'= 0.5e-7, 'interpolant'= false, 'output'= Array([eta]),  
            dsolveopts[]
         )
      )
   end proc,
   plotspec:= (Zremember, RX, RY),
   C:= plots:-contourplot(
      plotspec,
      #These default plot options can be changed by setting 'contouropts':
      'grid'= [25,25], 'contours'= 5, 'filled',
      'coloring'= ['yellow', 'orange'], 'color'= 'green',
      contouropts[]
   ),
   P:= plot3d(
      plotspec,
      #These default plot options can be changed by setting 'surfaceopts':
      'grid'= [25,25], 'style'= 'surfacecontour', 'contours'= 6,
      surfaceopts[]
   ),
   U, L #z-axis endpoints after margin adjustment
;
   #Stretch z-axis to include margins:
   (U,L):= ((Um,Lm,M,m)-> (M*(Lm-1)+m*Um, M*Lm+m*(Um-1)) /~ (Um+Lm-1))(
      zmargin[],
      (max,min)(op(3, indets(P, 'specfunc'('GRID'))[])) #actual z-axis range
   );
   plots:-display(  #final plot to display/return
      [
         plots:-spacecurve(
            {
               [[lhs(RX),rhs(RY),U],[rhs(RX),rhs(RY),U],[rhs(RX),rhs(RY),L]], #yz backwall
               [[rhs(RX),rhs(RY),U],[rhs(RX),lhs(RY),U],[rhs(RX),lhs(RY),L]]  #xz backwall
            },
            'color'= 'grey', 'thickness'= 0
         ),
         plottools:-transform((x,y)-> [x,y,L])(C), #dropped-shadow contours
         P
      ],
      #These default plot options can be changed simply by putting the option in the
      #ParamPlot3d call:
      'view'= ['DEFAULT', 'DEFAULT', L..U], 'orientation'= [-135, 75], 'axes'= 'frame',
      'labels'= [lhs(X), lhs(Y), Z], 'labelfont'= ['TIMES', 'BOLDOBLIQUE', 16],
      'caption'= nprintf(cat("%a = %4.2f, "$nops(FP)-1, "%a = %4.2f"), (lhs,rhs)~(FP)[]),
      'captionfont'= ['TIMES', 14],
      'projection'= 2/3,   
      _rest
   )
end proc:

ParamPlot3d(
   GetNu, S= 0..12, Ha= 0..12, [Rd= 4.5, Ec= 0.06, Nb= 0.15, Nt= 0.15, Sc= 0.5],
   labels= [S, Ha, Nu]
);

ParamPlot3d(
   GetNu, S= 0..12, Rd= 1..8, [Ha= 6, Ec= 0.06, Nb= 0.15, Nt= 0.15, Sc= 0.5],
   labels= [S, Rd, Nu]
);

ParamPlot3d(
   GetNu, S= 0..12, Ec= 0.01..0.10, [Ha= 6, Rd= 4.5, Nb= 0.15, Nt= 0.15, Sc= 0.5],
   labels= [S, Ec, Nu]
);

ParamPlot3d(
   GetNu, S= 0..12, Nt= 0.1..0.2, [Ha= 6, Rd= 4.5, Ec= 0.06, Nb= 0.15, Sc= 0.5],
   labels= [S, Nt, Nu]
);

ParamPlot3d(
   GetNu, Ha= 0..12, Rd= 1.0..8.0, [S= 6, Ec= 0.06, Nb= 0.15, Nt= 0.15, Sc= 0.5],
   labels= [Ha, Rd, Nu]
);

ParamPlot3d(
   GetNu, Ha= 0..12, Ec= 0.01..0.1, [S= 6, Rd= 4.5, Nb= 0.15, Nt= 0.15, Sc= 0.5],
   labels= [Ha, Ec, Nu]
);
   

ParamPlot3d(
   GetNu, Ha= 0..12, Nt= 0.1..0.2, [S= 6, Rd= 4.5, Ec= 0.06, Nb= 0.15, Sc= 0.5],
   labels= [Ha, Nt, Nu]
);

ParamPlot3d(
   GetNu, Ha= 0..12, Nb= 0.1..0.2, [S= 6, Rd= 4.5, Ec= 0.06, Nt= 0.15, Sc= 0.5],
   labels= [Ha, Nb, Nu]
);

ParamPlot3d(
   GetNu, Rd= 1.0..8.0, Nt= 0.1..0.2, [S= 6, Ha= 6, Ec= 0.06, Nb= 0.15, Sc= 0.5],
   labels= [Rd, Nt, Nu]
);

ParamPlot3d(
   GetNu, Rd= 1.0..8., Nb= 0.1..0.2, [S= 6, Ha= 6, Ec= 0.06, Nt= 0.15, Sc= 0.5],
   labels= [Rd, Nb, Nu]
);

ParamPlot3d(
   GetNu, Ec= 0.01..0.1, Nt= 0.1..0.2, [S= 6, Ha= 6, Rd= 4.5, Nb= 0.15, Sc= 0.5],
   labels= [Ec, Nt, Nu]
);

ParamPlot3d(
   GetNu, Ec= 0.01..0.1, Nb= 0.1..0.2, [S= 6, Ha= 6, Rd= 4.5, Nt= 0.15, Sc= 0.5],
   labels= [Ec, Nb, Nu]
);

#parameter values that remain fixed for the entire set of plots:
Pc:= Sc= 0.5, Ec= 0.1, Pr= 10:
 

#parameter values that remain fixed within each plot:
Ps:= [
   [S= 0.5, Ha= 2, Nt= 0.1, Nb= 0.1],
   [Ha= 2, Rd= 0.5, Nt= 0.1, Nb= 0.1],
   [S= 0.5, Rd= 0.5, Nt= 0.1, Nb= 0.1],
   [S= 0.5, Ha= 2, Rd= 0.5, Nb= 0.1],
   [S= 0.5, Ha= 2, Rd= 0.5, Nt= 0.1]
]:

#parameter value for each curve in each plot
Pv:= [
   Rd = [0.5, 2, 4, 8],
   S= [2, 4, 8, 12],
   Ha= [2, 4, 8, 12],
   Nt= [0.01, 0.05, 0.1, 0.2],
   Nb= [0.1, 0.2, 0.4, 0.6]
]:
      

for i to nops(Ps) do
   plots:-display(
      [seq(
         plots:-odeplot(
            Solve(lhs(Pv[i])= rhs(Pv[i])[j], Ps[i][], Pc),
            [eta, phi(eta)], 'color'= colseq[j], 'legend'= [lhs(Pv[i])= rhs(Pv[i])[j]]
         ),
         j= 1..nops(rhs(Pv[i]))
      )],
      'axes'= 'boxed', 'gridlines'= false,
      'labelfont'= ['TIMES', 'BOLDOBLIQUE', 16],
      'caption'= nprintf(
         cat("\n%a = %4.2f, "$nops(Ps[i])-1, "%a = %4.2f\n\n"), (lhs,rhs)~(Ps[i])[]
       ),
       'captionfont'= ['TIMES', 16]
   )
od;

NusseltList:= [entries(Solve:-SavedData, nolist)]:

nops(NusseltList);

V:= [S, Ha, Rd, Ec, Nt, Nb]:

NusseltData:= Matrix(map(R-> [map(v-> R[v], V)[], R[Nu]], NusseltList));