(************** Content-type: application/mathematica **************
                     CreatedBy='Mathematica 4.2'

                    Mathematica-Compatible Notebook

This notebook can be used with any Mathematica-compatible
application, such as Mathematica, MathReader or Publicon. The data
for the notebook starts with the line containing stars above.

To get the notebook into a Mathematica-compatible application, do
one of the following:

* Save the data starting with the line of stars above into a file
  with a name ending in .nb, then open the file inside the
  application;

* Copy the data starting with the line of stars above to the
  clipboard, then use the Paste menu command inside the application.

Data for notebooks contains only printable 7-bit ASCII and can be
sent directly in email or through ftp in text mode.  Newlines can be
CR, LF or CRLF (Unix, Macintosh or MS-DOS style).

NOTE: If you modify the data for this notebook not in a Mathematica-
compatible application, you must delete the line below containing
the word CacheID, otherwise Mathematica-compatible applications may
try to use invalid cache data.

For more information on notebooks and Mathematica-compatible 
applications, contact Wolfram Research:
  web: http://www.wolfram.com
  email: info@wolfram.com
  phone: +1-217-398-0700 (U.S.)

Notebook reader applications are available free of charge from 
Wolfram Research.
*******************************************************************)

(*CacheID: 232*)


(*NotebookFileLineBreakTest
NotebookFileLineBreakTest*)
(*NotebookOptionsPosition[     69023,       1713]*)
(*NotebookOutlinePosition[     70121,       1748]*)
(*  CellTagsIndexPosition[     70077,       1744]*)
(*WindowFrame->Normal*)



Notebook[{
Cell[TextData[StyleBox["Purely Incremental Solution Methods For Nonlinear \
Structural Analysis \nVersion: March 31, 2012 ",
  FontSize->14]], "Section",
  CellFrame->True,
  CellMargins->{{14, 37}, {Inherited, Inherited}},
  CellLabelMargins->{{8, Inherited}, {Inherited, Inherited}},
  ImageRegion->{{-0, 1}, {0, 1}},
  Background->RGBColor[1, 1, 0]],

Cell[TextData[{
  StyleBox["This Notebook is a testbed of incremental solution methods, \
converted to ",
    FontSize->12],
  StyleBox["Mathematica ",
    FontSize->12,
    FontSlant->"Italic"],
  StyleBox["from\nancient Fortran (1972-97).  Development history at Lockheed \
Palo Alto & CU Boulder:\n\n.  1966 PI Fotran 4 code (IBM 7094) for thesis \
(geometric + material nonlinearities) \n",
    FontSize->12],
  StyleBox[".", "Section",
    FontSize->12],
  StyleBox["  1972-73 Fortran code (Univac) for submerged antennas under \
steady ocean currents\n           Incorporated corrective phase and \
singularity regularization\n",
    FontSize->12],
  StyleBox[". ", "Section",
    FontSize->12],
  StyleBox[" 1980-82 Expanded Vax version in f77 for T2 Navy system - \
incorporated \n   Riks' arclength control. Code reached ~ 150000 lines\n",
    FontSize->12],
  StyleBox[". ", "Section",
    FontSize->12],
  StyleBox[" 1989  Tiny portion (about 1000 f77 lines) extracted to support \
NFEM first offering @ CU.\n   Removed element library, sparse solver, control \
language, database manager.\n",
    FontSize->12],
  StyleBox[". ", "Section",
    FontSize->12],
  StyleBox[" 1997  Converted to ",
    FontSize->12],
  StyleBox["Mathematica",
    FontSize->12,
    FontSlant->"Italic"],
  StyleBox[" 3 but split into benchmark-specific codes. Removed\n   \
correction phase, acceleration and line search.  Corrections restored in a \
subsequent Chapter.\n",
    FontSize->12],
  StyleBox[".", "Section",
    FontSize->12],
  StyleBox["  2012  (This notebook) Benchmarks gathered under single code - \
also unified \n   advancing methods & increment constraints, and provided \
better online display\n      \nMethods are applied to benchmark problems \
described in the newer Chapter 21 of the \nNFEM course.  Presently \
implemented benchmarks: MT and CG.\n\nPresent implementation include seven  \
RK-type step advancing integrators:\n    FE:                Forward Euler\n   \
 EMR:            ExplicitMidpoint Rule (a 2-stage Runge Kutta)\n    ETR:      \
       ExplicitTrapezoidal Rule  (a 2-stage Runge Kutta)\n    RK3R,RK3H:    \
3-stage Runge Kutta, two variants: R for Ralston, H for Heun\n    RK4C,RK4K:  \
  4-stage Runge Kutta, two variants: C for classic, K for Kutta\nWhich are \
combined with four Increment Control Strategies:\n    LC:         Lambda \
Control, also called Load Control\n    SC:         State Control, using a \
norm of the velocity vector\n    AC :        ArcLength Control  \
(Riks-Wempner)\n    HC :       Hyperelliptic Control (Felippa-Skeie)\n The \
implementations below use a constant stepsize ell and a positive-work \
criterion to\n traverse limit and turning points.  No  provision is made for \
bifurcation points.\n \n To run a problem driver script, all cells marked as \
\"initialization cells\"  ",
    FontSize->12],
  StyleBox["must be initialized",
    FontSize->12,
    FontSlant->"Italic"],
  StyleBox[". \n Those include all cells from the first one up to, and \
including, ",
    FontSize->12],
  StyleBox["NonLinSolver",
    FontSize->12,
    FontWeight->"Bold"]
}], "Text",
  CellFrame->True,
  CellMargins->{{14, 37}, {Inherited, Inherited}},
  CellLabelMargins->{{8, Inherited}, {Inherited, Inherited}},
  TextJustification->1,
  ImageRegion->{{-0, 1}, {0, 1}},
  FontSize->14,
  Background->RGBColor[1, 0.804959, 0.350027]],

Cell[TextData[StyleBox["Miscellaneous utilities",
  FontWeight->"Bold",
  FontColor->GrayLevel[0]]], "Text",
  CellFrame->True,
  CellMargins->{{14, 37}, {Inherited, Inherited}},
  CellLabelMargins->{{8, Inherited}, {Inherited, Inherited}},
  ImageRegion->{{-0, 1}, {0, 1}},
  Background->RGBColor[0, 1, 1]],

Cell[TextData[{
  StyleBox[" ",
    FontFamily->"Palatino"],
  "Mathematica version-dependent settings (to get plots displayed correctly)"
}], "Text",
  CellFrame->True,
  CellMargins->{{14, 37}, {Inherited, Inherited}},
  CellLabelMargins->{{8, Inherited}, {Inherited, Inherited}},
  ImageRegion->{{-0, 1}, {0, 1}},
  Background->RGBColor[1, 1, 0]],

Cell["\<\
Off[General::spell]; Off[General::spell1];
DisplayChannel=$DisplayFunction;
If [$VersionNumber>=6.0, DisplayChannel=Print]; 
  (* fix for Mathematica 6 & later *)\
\>", "Input",
  CellFrame->True,
  CellMargins->{{14, 37}, {Inherited, Inherited}},
  InitializationCell->True,
  CellLabelMargins->{{8, Inherited}, {Inherited, Inherited}},
  ImageRegion->{{-0, 1}, {0, 1}},
  Background->RGBColor[0, 1, 0]],

Cell[TextData[{
  StyleBox["ButcherTable:  ",
    FontWeight->"Bold",
    FontColor->GrayLevel[0]],
  StyleBox["returns specs of implemented RK-style integrators using the \
Butcher table.",
    FontColor->GrayLevel[0],
    FontVariations->{"CompatibilityType"->0}],
  StyleBox["\nFor details on what this table means, see any advanced book on \
numerical integration of ODE.",
    FontColor->GrayLevel[0]]
}], "Text",
  CellFrame->True,
  CellMargins->{{14, 37}, {Inherited, Inherited}},
  CellLabelMargins->{{8, Inherited}, {Inherited, Inherited}},
  ImageRegion->{{-0, 1}, {0, 1}},
  Background->RGBColor[1, 1, 0]],

Cell[TextData[{
  StyleBox["ButcherTable[int_]:=Module[{A,c,w},\n   If [int==\"FE\",\n       \
A={{0}}; c={0}; w={1}; \n       Return[{A,c,w}]]; \n   If [int==\"EMR\",\n    \
   A={{0,0},{1/2,0}}; c={0,1/2}; w={0,1}; \n       Return[{A,c,w}]];\n   If \
[int==\"ETR\",\n       A={{0,0},{1,0}}; c={0,1}; w={1/2,1/2}; \n       \
Return[{A,c,w}]]; \n   If [int==\"RK3R\",\n       \
A={{0,0,0},{1/2,0,0},{0,3/4,0}}; \n       c={0,1/2,3/4}; w={2/9,1/3,4/9}; \
Return[{A,c,w}]]; \n   If [int==\"RK3H\",\n       \
A={{0,0,0},{1/3,0,0},{0,2/3,0}}; c={0,1/3,2/3}; \n       w={1/4,0,3/4}; \
Return[{A,c,w}]]; \n   If [int==\"RK4C\",\n       \
A={{0,0,0,0},{1/2,0,0,0},{0,1/2,0,0},{0,0,1,0}}; \n       c={0,1/2,1/2,1}; \
w={1/6,1/3,1/3,1/6}; \n       Return[{A,c,w}]];\n   If [int==\"RK4K\",\n      \
 A={{0,0,0,0},{1/3,0,0,0},{-1/3,1,0,0},{1,-1,1,0}}; \n       c={0,1/3,2/3,1}; \
w={1/8,3/8,3/8,1/8}; \n       Return[{A,c,w}]];\n   Print[\"ButcherTable: \
Integrator \",int,\" not found\"];\n   Return[{Null,Null,Null}]];\n",
    FontColor->GrayLevel[0]],
  StyleBox["\n\
(*intab={\"FE\",\"EMR\",\"ETR\",\"RK3R\",\"RK3H\",\"RK4C\",\"RK4K\"};  \nFor \
[i=1,i<=Length[intab],i++, int=intab[[i]]; \n    {A,c,w}=ButcherTable[int];\n \
   Print [\"Integ: \",int,\" A=\",A//MatrixForm,\n     \" c=\",c//MatrixForm,\
\" w=\",w]];*)",
    FontColor->RGBColor[0, 0, 1]]
}], "Input",
  CellFrame->True,
  CellMargins->{{14, 37}, {Inherited, Inherited}},
  InitializationCell->True,
  CellLabelMargins->{{8, Inherited}, {Inherited, Inherited}},
  ImageRegion->{{-0, 1}, {0, 1}},
  Background->RGBColor[0, 1, 0]],

Cell[TextData[{
  StyleBox["MatrixFrobNorm",
    FontWeight->"Bold"],
  ": returns",
  StyleBox[" the Frobenius norm of a real square matrix",
    FontColor->GrayLevel[0]]
}], "Text",
  CellFrame->True,
  CellMargins->{{14, 37}, {Inherited, Inherited}},
  CellLabelMargins->{{8, Inherited}, {Inherited, Inherited}},
  ImageRegion->{{-0, 1}, {0, 1}},
  Background->RGBColor[1, 1, 0]],

Cell[TextData[{
  "MatrixFrobNorm[A_]:=Module[{i,j,n=Length[A],Fnorm=0},\n   If \
[n<=0,Return[0]];\n   For [i=1,i<=n,i++, Fnorm+=A[[i]].A[[i]]];\n   \
Return[Sqrt[Fnorm]/n]];",
  StyleBox["   ",
    FontColor->RGBColor[0, 0, 1]]
}], "Input",
  CellFrame->True,
  CellMargins->{{14, 37}, {Inherited, Inherited}},
  InitializationCell->True,
  CellLabelMargins->{{8, Inherited}, {Inherited, Inherited}},
  ImageRegion->{{-0, 1}, {0, 1}},
  Background->RGBColor[0, 1, 0]],

Cell[TextData[{
  StyleBox[" ICSFactor",
    FontWeight->"Bold"],
  ":  returns the factor",
  StyleBox[" f",
    FontSlant->"Italic"],
  " associated with the given incremental control strategy.\nArbitrary \
scaling of both control parameter \[Lambda] as well as velocity vector ",
  StyleBox["v",
    FontWeight->"Bold"],
  " entries is \npossible using the cscale (constraint scale) argument"
}], "Text",
  CellFrame->True,
  CellMargins->{{14, 37}, {Inherited, Inherited}},
  CellLabelMargins->{{8, Inherited}, {Inherited, Inherited}},
  ImageRegion->{{-0, 1}, {0, 1}},
  Background->RGBColor[1, 1, 0]],

Cell[TextData[{
  "ICSFactor[",
  StyleBox["ics_,",
    FontColor->GrayLevel[0]],
  "cscale",
  StyleBox["_,v_,q_]:=Module[{sgn=Sign[q.v],\n   \
csc=cscale,\[Lambda]scale=0,vscale=0,vv},\n   If [ics==\"LC\",Return[sgn]]; \n\
   If [Head[csc]==List&&Length[csc]<=0,csc=1];\n",
    FontColor->GrayLevel[0]],
  StyleBox[" (*Print[\"ics=\",ics,\" cscale=\",cscale,\n  \" \
len=\",Length[cscale],\" head=\",Head[cscale]];*)",
    FontColor->RGBColor[1, 0, 0]],
  StyleBox["\n   If [Length[",
    FontColor->GrayLevel[0]],
  "csc",
  StyleBox["]<=0,vscale=",
    FontColor->GrayLevel[0]],
  "csc",
  StyleBox["];\n   If [Length[",
    FontColor->GrayLevel[0]],
  "csc",
  StyleBox["]==1,vscale=",
    FontColor->GrayLevel[0]],
  "csc",
  StyleBox["[[1]]];\n   If [Length[",
    FontColor->GrayLevel[0]],
  "csc",
  StyleBox["]==2,{\[Lambda]scale,vscale}=",
    FontColor->GrayLevel[0]],
  "csc",
  StyleBox["];\n   If [Head[vscale]==List, vv=(vscale.v)^2,\n       \
vv=(vscale*v).(vscale*v),vv=(vscale*v).(vscale*v)];\n   If \
[ics==\"SC\",Return[sgn*Sqrt[vv]]];  \n   If \
[ics==\"AC\",Return[sgn*Sqrt[1+vv]]];\n   If [ics==\"HC\",Return[sgn*Sqrt[\
\[Lambda]scale^2+vv]]];\n   Return[Null]];",
    FontColor->GrayLevel[0]],
  StyleBox[" \n\n(*v={v1,v2,v3}; q={q1,q2,q3}; \n\
Print[ICSFactor[\"LC\",{},v,q]];  \nPrint[ICSFactor[\"SC\",{},v,q]]; \n\
Print[ICSFactor[\"SC\",1,v,q]]; \nPrint[ICSFactor[\"SC\",{{0,1,0}},v,q]]; \n\
Print[ICSFactor[\"SC\",{{2,3,1/4}},v,q]]; \nPrint[ICSFactor[\"AC\",{},v,q]]; \
\nPrint[ICSFactor[\"AC\",1,v,q]]; \nPrint[ICSFactor[\"AC\",{{0,1,0}},v,q]]; \n\
Print[ICSFactor[\"AC\",{{2,3,1/4}},v,q]]; \nPrint[ICSFactor[\"HC\",{},v,q]];\n\
Print[ICSFactor[\"HC\",{1,1},v,q]]; \n\
Print[ICSFactor[\"HC\",{1/2,{0,1,0}},v,q]]; \n\
Print[ICSFactor[\"HC\",{6,{2,3,1/4}},v,q]]; *)\n",
    FontColor->RGBColor[0, 0, 1]]
}], "Input",
  CellFrame->True,
  CellMargins->{{14, 37}, {Inherited, Inherited}},
  InitializationCell->True,
  CellLabelMargins->{{8, Inherited}, {Inherited, Inherited}},
  ImageRegion->{{-0, 1}, {0, 1}},
  Background->RGBColor[0, 1, 0]],

Cell[TextData[{
  StyleBox["Display modules ",
    FontWeight->"Bold"],
  " support printing & plotting by driver scripts"
}], "Text",
  CellFrame->True,
  CellMargins->{{14, 37}, {Inherited, Inherited}},
  CellLabelMargins->{{8, Inherited}, {Inherited, Inherited}},
  ImageRegion->{{-0, 1}, {0, 1}},
  Background->RGBColor[0, 1, 1]],

Cell[TextData[{
  StyleBox["PrintPISolTabl",
    FontWeight->"Bold"],
  "e:  Module to print computed solution table.  Converted  from old \n\
Fortran 66 code (written in 1972 at Lockheed PARL) to ",
  StyleBox["Mathematica",
    FontSlant->"Italic"],
  " 4.1 in 1997.\nSignificant upgrade 2012. Nicer typesetting and more \
flexible (output cell may be \nmoved directly to documents via EPS) but \
limited to 6 DOF max. "
}], "Text",
  CellFrame->True,
  CellMargins->{{14, 37}, {Inherited, Inherited}},
  CellLabelMargins->{{8, Inherited}, {Inherited, Inherited}},
  ImageRegion->{{-0, 1}, {0, 1}},
  Background->RGBColor[1, 1, 0]],

Cell[TextData[{
  "PrintPISolTable[ptags_,soltab_,doflab_,digits_]:= \n  \
Module[{nsteps=Length[soltab],numdof,i,j,k,n,\[Lambda]n,un,vn,qn,\n    \
rn,Kdetn,elln,an,\[Kappa]n,\[Kappa]0,rem,vnorm,rnorm,t,th,thbeg,thend,\n    \
thstate,d\[Lambda],f\[Lambda],du,fu,dv,fv,dr,fr,dK,fK,da,fa,d\[Kappa],f\
\[Kappa],\n    dftab=Table[{5,3},{7}],diglen,choptol=10.0^(-12),\n    \
modname=\"PrintPISolTable]\",problem,int,ics,\n    \
TF,TFB,FF=\"Courier\",FS=10,FW=\"Plain\",FSL=\"Plain\"},\n  \
TF[text_]:=StyleForm[text,FontFamily->FF,FontSize->FS,\n            \
FontWeight->FW,FontSlant->FSL];\n  \
TFB[text_]:=StyleForm[text,FontFamily->FF,FontSize->FS,\n            \
FontWeight->\"Bold\",FontSlant->FSL];\n  \
TFI[text_]:=StyleForm[text,FontFamily->\"Times\",FontSize->FS,\n            \
FontWeight->FW,FontSlant->\"Italic\"];\n  If [nsteps<=0, Print [modname,\": \
soltab empty. Print skipped.\"]; \n      Return[]]; \
numdof=Length[soltab[[1,3]]]; \n  If [numdof<=0||numdof>6, Print [modname,\": \
only 1-6 DOF allowed.\",\n     \" Print skipped.\"]; Return[]];\n  t=Table[\" \
\",{nsteps},{numdof+8}]; diglen=Length[digits];\n  If [diglen==1, \
dftab=Table[digits[[1]],{7}]];\n  If [diglen==7, dftab=digits];\n  {{d\
\[Lambda],f\[Lambda]},{du,fu},{dv,fv},{dr,fr},{dK,fK},{da,fa},{d\[Kappa],f\
\[Kappa]}}=dftab; \n  For [i=1,i<=nsteps,i++,\n     \
{n,\[Lambda]n,un,vn,rn,qn,Kdetn,elln,an,\[Kappa]n,rem}=soltab[[i]];\n     \
vnorm=Sqrt[vn.vn]; rnorm=Sqrt[rn.rn]; \n     \
{\[Lambda]n,un,vnorm,rnorm,Kdetn,elln,an,\[Kappa]n}=  Chop[N[\n     \
{\[Lambda]n,un,vnorm,rnorm,Kdetn,elln,an,\[Kappa]n}],choptol];\n        \
t[[i,1]]=TF[ToString[n]];\n        t[[i,2]]=TFB[PaddedForm[\[Lambda]n,{d\
\[Lambda],f\[Lambda]}]];\n        For [j=1,j<=numdof,j++,\n             \
t[[i,j+2]]=TFB[PaddedForm[un[[j]],{du,fu}]]]; \n        \
t[[i,numdof+3]]=TF[PaddedForm[vnorm,{dv,fv}]];\n        \
t[[i,numdof+4]]=TF[PaddedForm[rnorm,{dr,fr}]];\n        \
t[[i,numdof+5]]=TF[PaddedForm[Kdetn,{dK,fK}]];\n        \
t[[i,numdof+6]]=TF[PaddedForm[an,{da,fa}]];\n        \
t[[i,numdof+7]]=TF[PaddedForm[\[Kappa]n,{d\[Kappa],f\[Kappa]}]];\n        \
t[[i,numdof+8]]=TFI[rem]]; \n   thend={\"v-norm\",\"r-norm\",\"Kdet \",\"a \
\",\"\[Kappa] \",\"Remarks\"};\n   thstate={\"u1  \",\"u2  \",\"u3  \",\"u4  \
\",\"u5  \",\"u6  \"};\n   thbeg={\"n\",\"\[Lambda] \"}; k=Length[doflab];\n  \
 If [k>0, For[i=1,i<=k,i++,thstate[[i]]=doflab[[i]] ]];\n   \
th=TF/@Flatten[Join[thbeg,Take[thstate,numdof],thend]];",
  StyleBox[" \n",
    FontColor->RGBColor[1, 0, 0]],
  "   {problem,int,ics}=ptags;\n   Print[TF[\"PI solution for problem: \
\"],TF[problem],\n          TF[\", with int: \"],TF[int],TF[\" & ics: \
\"],TF[ics]]; \n   Print[TableForm[t,TableAlignments->{Right},\n         \
TableDirections->{Column,Row},TableSpacing->{0,1}, \n         TableHeadings \
->{None,th}]];\n   ClearAll[t];\n   ];  \n",
  StyleBox["\nsoltab=\n{ {0,0.0,{ \
0.0,0},{0.02,0.05},{0.42,0},{0,-1},45.,0.02,1.0,1.0,\"ref config\"}, \n  \
{1,0.2,{-0.3,0},{0.03,-.04},{0.35,0},{0,-1},43.,0.02,1.0,0.9,\"step 1\"} };\n\
PrintPISolTable[{\"MT\",\"FE\",\"LC\"},soltab,{\"uX \",\"uY \"},{4,2}];\n\
soltab=\n{ {0,0.0,{ 0.0},{0.07},{0.42},{1},6.*10^8,0.02,1.0,1.0,\"ref \
config\"}, \n  {1,0.2,{-0.3},{-.12},{0.35},{1},5.*10^8,0.02,1.0,0.9,\"step \
1\"} };\nPrintPISolTable[{\"CG\",\"FE\",\"AC\"},soltab,{\"\[Mu]  \"}, {}];",
    FontColor->RGBColor[0, 0, 1]]
}], "Input",
  CellFrame->True,
  CellMargins->{{14, 37}, {Inherited, Inherited}},
  InitializationCell->True,
  CellLabelMargins->{{8, Inherited}, {Inherited, Inherited}},
  ImageRegion->{{-0, 1}, {0, 1}},
  Background->RGBColor[0, 1, 0]],

Cell[TextData[{
  StyleBox["PlotPISolTabl",
    FontWeight->"Bold"],
  "e:  Module to plot selected histories from solution table"
}], "Text",
  CellFrame->True,
  CellMargins->{{14, 37}, {Inherited, Inherited}},
  CellLabelMargins->{{8, Inherited}, {Inherited, Inherited}},
  ImageRegion->{{-0, 1}, {0, 1}},
  Background->RGBColor[1, 1, 0]],

Cell[TextData[{
  "PlotPISolTable[ptags_,soltab_,xytags_,doflab_,imgsiz_]:= \n  \
Module[{vartag,varidx,nsol,xvar,yvar,ntag,\n   \
i,ix,iy,j,k,m,n,ixlen,iylen,t,xyv,title,\n   \
problem,int,ics,modname=\"PlotPISolTable:\"},\n  vartag={\"step\",\"n\",\"\
\[Lambda]\",\"u1\", \"u2\", \"u3\", \"u4\", \"u5\", \"u6\",\n  \
\"||u||\",\"||v||\",\"||r||\",\"||q||\",\"Kdet\",\"ell\",\"a\",\"\[Kappa]\"};\
\n  varidx={   1,   1,  2,{3,1},{3,2},{3,3},{3,4},{3,5},{3,6}, \n    -3,     \
-4,     -5,     -6,      7,    8,   9, 10};\n {xvar,yvar}=xytags; \
k=Length[doflab]; \n  If [k>0, \
For[i=1,i<=k,i++,vartag[[i+3]]=doflab[[i]]]];",
  StyleBox[" ",
    FontColor->RGBColor[1, 0, 0]],
  "\n  ntag=Length[vartag]; ix=iy=0;\n  For [i=1,i<=ntag,i++,\n       If \
[xvar==vartag[[i]], ix=varidx[[i]];Break[]]];\n  If [ix==0, Print[modname,\" \
unknown x-var tag \",xvar]; \n       Return[]];\n  For [i=1,i<=ntag,i++, \n   \
    If [yvar==vartag[[i]], iy=varidx[[i]];Break[]]];\n  If [iy==0, \
Print[modname,\" unknown y-var tag \",yvar]; \n       Return[]]; \n  \
ixlen=Length[ix]; iylen=Length[iy]; nsol=Length[soltab]; \n  \
m=Length[soltab[[1]]]; xyv=Table[{0.,0.},{nsol}];\n  For [n=1,n<=nsol,n++, \
t=N[Take[soltab[[n]],m-1]];\n       If [ixlen==0, k=Abs[ix];\n           If \
[ix>0, xyv[[n,1]]=t[[ix]]];\n           If [ix<0, \
xyv[[n,1]]=Sqrt[t[[k]].t[[k]]] ]];\n       If [ixlen==2, {i,j}=ix; \
xyv[[n,1]]=t[[i,j]]];\n       If [iylen==0, k=Abs[iy];\n           If [iy>0, \
xyv[[n,2]]=t[[iy]]];\n           If [iy<0, xyv[[n,2]]=Sqrt[t[[k]].t[[k]]] ]];\
\n       If [iylen==2, {i,j}=iy; xyv[[n,2]]=t[[i,j]]];\n       ]; \n  \
{problem,int,ics}=ptags;\n  title=yvar<>\" vs \"<>xvar<>\" for problem \
\"<>problem<>\n        \", int: \"<>int<>\", ics: \"<>ics;\n  \
ListPlot[xyv,PlotJoined->True,PlotRange->All,\n    \
Frame->True,AxesOrigin->{0,0},ImageSize->imgsiz, \n    \
PlotStyle->{AbsoluteThickness[1.5],RGBColor[1,0,0]},\n    \
PlotLabel->title,DisplayFunction->DisplayChannel];\n    ClearAll[xyv];\n   ]; \
 \n\n(*",
  StyleBox["ptags={\"MT\",\"EMR\",\"AC\"}; imgsiz=300;\nsoltab=\n{{0,0,    \
{0,0},      {0.,-1.414},{0.,0.},     {0,-1.},0.5,  0.04,1,1,\" \"},\n \
{1,0.023,{0.,-0.032},{0.,-1.414},{0.,0.00112},{0,-1.},0.5,  0.04,1,1,\" \"},\n\
 {2,0.045,{0.,-0.066},{0.,-1.565},{0.,0.00227},{0,-1.},0.437,0.04,1,1,\" \"},\
\n {3,0.064,{0.,-0.101},{0.,-1.751},{0.,0.00345},{0,-1.},0.378,0.04,1,1,\" \
\"}};\n PrintPISolTable[ptags,soltab,",
    FontColor->RGBColor[0, 0, 1]],
  StyleBox["{\"uX\",\"uY\"}",
    FontColor->RGBColor[0, 0, 1]],
  StyleBox[",{}];\n PlotPISolTable [ptags,soltab,{\"uY\",\"\[Lambda]\"},",
    FontColor->RGBColor[0, 0, 1]],
  StyleBox["{\"uX\",\"uY\"},imgsiz",
    FontColor->RGBColor[0, 0, 1]],
  StyleBox["];\n PlotPISolTable [ptags,soltab,{\"||r||\",\"\[Lambda]\"},",
    FontColor->RGBColor[0, 0, 1]],
  StyleBox["{\"uX\",\"uY\"},imgsiz",
    FontColor->RGBColor[0, 0, 1]],
  StyleBox["];*)",
    FontColor->RGBColor[0, 0, 1]]
}], "Input",
  CellFrame->True,
  CellMargins->{{14, 37}, {Inherited, Inherited}},
  InitializationCell->True,
  CellLabelMargins->{{8, Inherited}, {Inherited, Inherited}},
  ImageRegion->{{-0, 1}, {0, 1}},
  Background->RGBColor[0, 1, 0]],

Cell[TextData[{
  StyleBox["Benchmark problem specific  modules.  ",
    FontWeight->"Bold"],
  StyleBox["Two implemented so far: MT and CG",
    FontVariations->{"CompatibilityType"->0}]
}], "Text",
  CellFrame->True,
  CellMargins->{{14, 37}, {Inherited, Inherited}},
  CellLabelMargins->{{8, Inherited}, {Inherited, Inherited}},
  ImageRegion->{{-0, 1}, {0, 1}},
  Background->RGBColor[0, 1, 1]],

Cell[TextData[{
  StyleBox["TotalResidualOfMisesTruss, IntForceOfMisesTruss, \
ExtForceOfMisesTruss,\nTanStiffOfMisesTruss, IncLoadOfMisesTruss,  \
TanStiffDetOfMisesTruss:",
    FontWeight->"Bold"],
  "\ntotal residual force vector, internal force vector, external force \
vector,\ntangent stiffness matrix, incremental load vector and tangent \
stiffness\ndeterminant, respectively, of von Mises Truss benchmark"
}], "Text",
  CellFrame->True,
  CellMargins->{{14, 37}, {Inherited, Inherited}},
  CellLabelMargins->{{8, Inherited}, {Inherited, Inherited}},
  ImageRegion->{{-0, 1}, {0, 1}},
  Background->RGBColor[1, 1, 0]],

Cell[TextData[{
  "TotalResidualOfMisesTruss[MTprop_,\[Lambda]_,state_,numer_]:=\n   \
Module[{pv,fv,rv}, \n   \
pv=IntForceOfMisesTruss[MTprop,\[Lambda],state,numer]; \n   \
fv=ExtForceOfMisesTruss[MTprop,\[Lambda],state,numer]; \n   rv=pv-fv; If \
[!numer, rv=Simplify[rv]];\n   Return[rv]];\n\nIntForceOfMisesTruss[MTprop_,\
\[Lambda]_,state_,numer_]:=\n   Module[{S,H,Em,A0,kP,uX,uY,c,p}, \n   If \
[!numer, {S,H,Em,A0,kP}=MTprop;    {uX,uY}=state]; \n   If [ numer, \
{S,H,Em,A0,kP}=N[MTprop]; {uX,uY}=N[state]]; \n   \
c=4*Em*A0/(4*H^2+S^2)^(3/2);\n   p=c*{uX*(S^2+2*uX^2+4*H*uY+2*uY^2),\n        \
2*(H+uY)*(uX^2+2*H*uY+uY^2)};\n   Return[p]];\n     \n\
ExtForceOfMisesTruss[MTprop_,\[Lambda]_,state_,numer_]:=\n   Module[{f}, f=\
\[Lambda]*IncLoadOfMisesTruss[MTprop,\[Lambda],state,numer];\n   If [!numer, \
f=Simplify[f]]; \n   Return[f]];\n\n\
TanStiffOfMisesTruss[MTprop_,\[Lambda]_,state_,numer_]:=\n   \
Module[{S,H,Em,A0,kP,uX,uY,c,K}, \n   If [!numer, {S,H,Em,A0,kP}=MTprop;    \
{uX,uY}=state]; \n   If [ numer, {S,H,Em,A0,kP}=N[MTprop]; {uX,uY}=N[state]]; \
\n   c=4*Em*A0/(4*H^2+S^2)^(3/2);\n   \
K=c*{{(S^2+6*uX^2+4*H*uY+2*uY^2),4*uX*(H+uY)},\n         \
{4*uX*(H+uY),2*(2*H^2+uX^2+6*H*uY+3*uY^2)}};\n   Return[K]];\n  \n\
IncLoadOfMisesTruss[MTprop_,\[Lambda]_,state_,numer_]:=\n   \
Module[{S,H,Em,A0,kP,q}, \n   If [!numer, {S,H,Em,A0,kP}=MTprop]; \n   If [ \
numer, {S,H,Em,A0,kP}=N[MTprop]]; \n   q={0,-Em*A0};\n   Return[q]];\n     \n\
TanStiffDetOfMisesTruss[MTprop_,\[Lambda]_,state_,numer_]:=\n   \
Module[{K,Kdet}, K=TanStiffOfMisesTruss[MTprop,\[Lambda],state,numer];\n   \
Kdet=Det[K]; If [!numer, Kdet=Simplify[Kdet]];\n   Return[Kdet]];\n  ",
  StyleBox["   \n(*ClearAll[S,H,Em,A0,\[Lambda],uX,uY];  \
MTprop={S,H,Em,A0,0};\np= \
IntForceOfMisesTruss[MTprop,\[Lambda],{uX,uY},False];\nK= \
TanStiffOfMisesTruss[MTprop,\[Lambda],{uX,uY},False];\nKK={D[p,uX],D[p,uY]};\n\
Print[\"Check K=grad(p): \",Simplify[K-KK]];\n\
Kdet=TanStiffDetOfMisesTruss[MTprop,\[Lambda],{0,uY},False];\n\
Print[\"Kdet=\",Kdet];\n\
r=TotalResidualOfMisesTruss[MTprop,\[Lambda],{uX,uY},False];\n\
Print[\"r=\",r//MatrixForm];*)",
    FontColor->RGBColor[0, 0, 1]]
}], "Input",
  CellFrame->True,
  CellMargins->{{14, 37}, {Inherited, Inherited}},
  InitializationCell->True,
  CellLabelMargins->{{8, Inherited}, {Inherited, Inherited}},
  ImageRegion->{{-0, 1}, {0, 1}},
  Background->RGBColor[0, 1, 0]],

Cell[TextData[{
  StyleBox["TotalResidualOfCircleGame,  TanStiffOfCircleGame, \
IncLoadOfCircleGame,  \nTanStiffDetOfCircleGame: ",
    FontWeight->"Bold"],
  "total residual force vector, tangent stiffness matrix, \nincremental load \
vector and tangent stiffness determinant, respectively, of Circle Game \n\
benchmark (named in honor of Joni Mitchell's '66 \"carousel of time\" song - \
\ntype \"Joni Mitchell Circle Game\" on YouTube)"
}], "Text",
  CellFrame->True,
  CellMargins->{{14, 37}, {Inherited, Inherited}},
  CellLabelMargins->{{8, Inherited}, {Inherited, Inherited}},
  ImageRegion->{{-0, 1}, {0, 1}},
  Background->RGBColor[1, 1, 0]],

Cell[TextData[{
  "TotalResidualOfCircleGame[CGprop_,\[Lambda]_,state_,numer_]:=\n   \
Module[{bipath,\[Mu],r}, bipath=CGprop[[1]]; \[Mu]=state[[1]];\n   r={\
\[Lambda]^2+\[Mu]^2-1}; If [bipath, r=(\[Lambda]-\[Mu])*r];\n   If [!numer, \
r=Simplify[r]];\n   Return[r]]\n\n\
TanStiffOfCircleGame[CGprop_,\[Lambda]_,state_,numer_]:=\n   Module[{bipath,\
\[Mu],K}, bipath=CGprop[[1]]; \[Mu]=state[[1]];\n   K={{\[InvisibleSpace]2*\
\[Mu]}}; If [bipath, K={{\[InvisibleSpace]1-\[Lambda]^2+2*\[Lambda]*\[Mu]-3*\
\[Mu]^2}}];\n   If [!numer, K=Simplify[K]];\n   Return[K]];\n  \n\
IncLoadOfCircleGame[CGprop_,\[Lambda]_,state_,numer_]:=\n   Module[{bipath,\
\[Mu],q}, bipath=CGprop[[1]]; \[Mu]=state[[1]];\n   q=\[InvisibleSpace]-{2*\
\[Lambda]}; If [bipath, \
q={\[InvisibleSpace]1-3*\[Lambda]^2+2*\[Lambda]*\[Mu]-\[Mu]^2}];\n   If \
[!numer, q=Simplify[q]];\n   Return[q]];\n     \n\
TanStiffDetOfCircleGame[CGprop_,\[Lambda]_,state_,numer_]:=\n   \
Module[{K,Kdet}, \n   K=TanStiffOfCircleGame[CGprop,\[Lambda],state,numer];\n \
  Kdet=K; If [!numer, Kdet=Simplify[Kdet]];\n   Return[",
  "Kdet",
  "]];\n\n(*",
  StyleBox["ClearAll[\[Lambda],\[Mu],bipath]; bipath=True;   \nr= \
TotalResidualOfCircleGame[{bipath},\[Lambda],{\[Mu]},False]; Print[\"r=\",r];\
\nK= TanStiffOfCircleGame     [{bipath},\[Lambda],{\[Mu]},False]; \
Print[\"K=\",K//InputForm];\nq= IncLoadOfCircleGame      \
[{bipath},\[Lambda],{\[Mu]},False]; Print[\"q=\",q//InputForm];*)",
    FontColor->RGBColor[0, 0, 1]]
}], "Input",
  CellFrame->True,
  CellMargins->{{14, 37}, {Inherited, Inherited}},
  InitializationCell->True,
  CellLabelMargins->{{8, Inherited}, {Inherited, Inherited}},
  ImageRegion->{{-0, 1}, {0, 1}},
  Background->RGBColor[0, 1, 0]],

Cell[TextData[{
  StyleBox["ProblemInfo",
    FontWeight->"Bold"],
  ": benchmark problem data used by driver script and/or solver"
}], "Text",
  CellFrame->True,
  CellMargins->{{14, 37}, {Inherited, Inherited}},
  CellLabelMargins->{{8, Inherited}, {Inherited, Inherited}},
  ImageRegion->{{-0, 1}, {0, 1}},
  Background->RGBColor[1, 1, 0]],

Cell["\<\
ProblemInfo[problem_]:=Module[{numdof=0},
   If [problem==\"MT\", numdof=2];
   If [problem==\"CG\", numdof=1];
   Return[numdof]];\
\>", "Input",
  CellFrame->True,
  CellMargins->{{14, 37}, {Inherited, Inherited}},
  InitializationCell->True,
  CellLabelMargins->{{8, Inherited}, {Inherited, Inherited}},
  ImageRegion->{{-0, 1}, {0, 1}},
  Background->RGBColor[0, 1, 0]],

Cell[TextData[{
  StyleBox["Interface modules",
    FontWeight->"Bold"],
  StyleBox[" link nonlinear solution modules to benchmark problems",
    FontVariations->{"CompatibilityType"->0}]
}], "Text",
  CellFrame->True,
  CellMargins->{{14, 37}, {Inherited, Inherited}},
  CellLabelMargins->{{8, Inherited}, {Inherited, Inherited}},
  ImageRegion->{{-0, 1}, {0, 1}},
  Background->RGBColor[0, 1, 1]],

Cell[TextData[{
  StyleBox["TotalResidual",
    FontWeight->"Bold"],
  ": Benchmark interface module to get total residual vector "
}], "Text",
  CellFrame->True,
  CellMargins->{{14, 37}, {Inherited, Inherited}},
  CellLabelMargins->{{8, Inherited}, {Inherited, Inherited}},
  ImageRegion->{{-0, 1}, {0, 1}},
  Background->RGBColor[1, 1, 0]],

Cell[TextData[{
  "TotalResidual[ptags_,geopar_,matpar_,fabpar_,solpar_,\n         \
reffor_,conpar_,state_,options_]:=Module[{\n    S,H,Em,A0,A10,A20,A30,uX,uY,\
\[Lambda],\[Mu],numer,bipath,\n    problem,int,ics,pnlen,pnam2,pnam3,r},\n    \
{problem,int,ics}=ptags; pnlen=StringLength[problem];\n    \
pnam2=pnam3=StringTake[problem,{1,2}]; \n    If [pnlen>=3, \
pnam3=StringTake[problem,{1,3}]];\n    numer=options[[1]]; \
\[Lambda]=conpar[[1]];\n    If [pnam2==\"MT\", \n        {S,H}=geopar; \
{Em}=matpar; {A10,A20,A30}=fabpar; \n        A0=(A10+A20)/2;  {uX,uY}=state; \
\n        r=TotalResidualOfMisesTruss[{S,H,Em,A0,0},\[Lambda],{uX,uY},numer];\
\n        If [!numer,r=Simplify[r]]; Return[r]];\n    If [pnam2==\"CG\", \
\[Mu]=state[[1]]; bipath=False;\n        If [pnam3==\"CGB\", bipath=True]; \n \
       r=TotalResidualOfCircleGame[{bipath},\[Lambda],{\[Mu]},numer];\n       \
 If [!numer,r=Simplify[r]]; Return[r]];\n    Print[\" Problem \",problem,\" \
not implemented\"];\n    Return[Null]]; \n",
  StyleBox["    \n\
(*ClearAll[S,H,Em,A0,uX,uY,\[Lambda],\[Mu],numer,ptags,geopar,\n\
matpar,fabpar,solpar,",
    FontColor->RGBColor[0, 0, 1]],
  StyleBox["reffor",
    FontColor->RGBColor[0, 0, 1]],
  StyleBox[",conpar,state,options]; \nptags={\"MT\",\"FE\",\"LC\"}; \
geopar={S,H}; matpar={Em}; fabpar={A0,A0,0};  \nsolpar={}; \
conpar={\[Lambda]}; state={uX,uY}; numer=False; options={numer};\n\
r=TotalResidual[ptags,geopar,matpar,fabpar,solpar,\n  ",
    FontColor->RGBColor[0, 0, 1]],
  StyleBox["reffor",
    FontColor->RGBColor[0, 0, 1]],
  StyleBox[",conpar,state,options]; Print[\"r=\",r//MatrixForm];\n\
ptags={\"CGB\",\"FE\",\"LC\"}; state={\[Mu]}; \n\
r=TotalResidual[ptags,geopar,matpar,fabpar,solpar,\n  ",
    FontColor->RGBColor[0, 0, 1]],
  StyleBox["reffor",
    FontColor->RGBColor[0, 0, 1]],
  StyleBox[",conpar,state,options]; Print[\"r=\",r//MatrixForm]; *)",
    FontColor->RGBColor[0, 0, 1]]
}], "Input",
  CellFrame->True,
  CellMargins->{{14, 37}, {Inherited, Inherited}},
  InitializationCell->True,
  CellLabelMargins->{{8, Inherited}, {Inherited, Inherited}},
  ImageRegion->{{-0, 1}, {0, 1}},
  Background->RGBColor[0, 1, 0]],

Cell[TextData[{
  StyleBox["TanStiff",
    FontWeight->"Bold"],
  ": Benchmark interface module for tangent stiffness "
}], "Text",
  CellFrame->True,
  CellMargins->{{14, 37}, {Inherited, Inherited}},
  CellLabelMargins->{{8, Inherited}, {Inherited, Inherited}},
  ImageRegion->{{-0, 1}, {0, 1}},
  Background->RGBColor[1, 1, 0]],

Cell[TextData[{
  "TanStiff[ptags_,geopar_,matpar_,fabpar_,solpar_,\n         \
reffor_,conpar_,state_,options_]:=Module[{\n    S,H,Em,A0,A10,A20,A30,uX,uY,\
\[Lambda],\[Mu],numer,bipath,\n    problem,int,ics,pnlen,pnam2,pnam3,K},\n   \
{problem,int,ics}=ptags; pnlen=StringLength[problem];\n    \
pnam2=pnam3=StringTake[problem,{1,2}]; \n    If [pnlen>=3, \
pnam3=StringTake[problem,{1,3}]];\n    numer=options[[1]]; \
\[Lambda]=conpar[[1]];\n    If [pnam2==\"MT\", \n        {S,H}=geopar; \
{Em}=matpar; {A10,A20,A30}=fabpar; \n  ",
  StyleBox["      (*Print[\"S,H,Em,A10,A20,A30=\",{S,H,Em,A10,A20,A30}];*)",
    FontColor->RGBColor[1, 0, 0]],
  "\n        A0=(A10+A20)/2; {uX,uY}=state; \n      ",
  StyleBox["  (*Print[\" \[Lambda]=\",\[Lambda],\" uX,uY=\",{uX,uY}];*)",
    FontColor->RGBColor[1, 0, 0]],
  "\n        K=TanStiffOfMisesTruss[{S,H,Em,A0,0},\[Lambda],{uX,uY},numer];\n \
       If [!numer,K=Simplify[K]]; Return[K]];\n    If [pnam2==\"CG\", \
\[Mu]=state[[1]]; bipath=False;\n        If [pnam3==\"CGB\", bipath=True]; \n \
       K=TanStiffOfCircleGame[{bipath},\[Lambda],{\[Mu]},numer];\n        If \
[!numer,K=Simplify[K]]; Return[K]];\n    Print[\" Problem \",problem,\" not \
implemented\"];\n    Return[Null]];\n",
  StyleBox["    \n\
(*ClearAll[S,H,Em,A0,uX,uY,\[Lambda],\[Mu],numer,ptags,geopar,\n\
matpar,fabpar,solpar,reffor,conpar,state,options]; \n\
ptags={\"MT\",\"FE\",\"LC\"}; geopar={S,H}; matpar={Em}; fabpar={A0,A0,0};  \n\
solpar={}; conpar={\[Lambda]}; state={uX,uY}; numer=False; options={numer};\n\
r=TotalResidual[ptags,geopar,matpar,fabpar,solpar,\n      \
reffor,conpar,state,options]; Print[\"r=\",r//MatrixForm];  \n\
K=TanStiff[ptags,geopar,matpar,fabpar,solpar,\n      \
reffor,conpar,state,options]; Print[\"K=\",K//MatrixForm];\n\
ptags={\"CGB\",\"FE\",\"LC\"}; state={\[Mu]};  \n\
K=TanStiff[ptags,geopar,matpar,fabpar,solpar,\n      \
reffor,conpar,state,options]; Print[\"K=\",K//MatrixForm];*)",
    FontColor->RGBColor[0, 0, 1]]
}], "Input",
  CellFrame->True,
  CellMargins->{{14, 37}, {Inherited, Inherited}},
  InitializationCell->True,
  CellLabelMargins->{{8, Inherited}, {Inherited, Inherited}},
  ImageRegion->{{-0, 1}, {0, 1}},
  Background->RGBColor[0, 1, 0]],

Cell[TextData[{
  StyleBox["IncLoad",
    FontWeight->"Bold"],
  ": Benchmark interface module for incremental load vector "
}], "Text",
  CellFrame->True,
  CellMargins->{{14, 37}, {Inherited, Inherited}},
  CellLabelMargins->{{8, Inherited}, {Inherited, Inherited}},
  ImageRegion->{{-0, 1}, {0, 1}},
  Background->RGBColor[1, 1, 0]],

Cell[TextData[{
  "IncLoad[ptags_,geopar_,matpar_,fabpar_,solpar_,\n        \
reffor_,conpar_,state_,options_]:= Module[{\n    S,H,Em,A0,A10,A20,A30,uX,uY,\
\[Lambda],\[Mu],numer,bipath,\n    problem,int,ics,pnlen,pnam2,pnam3,q},\n   \
{problem,int,ics}=ptags; pnlen=StringLength[problem];\n    \
pnam2=pnam3=StringTake[problem,{1,2}]; \n    If [pnlen>=3, \
pnam3=StringTake[problem,{1,3}]];\n    numer=options[[1]]; \
\[Lambda]=conpar[[1]];\n    If [pnam2==\"MT\", \n        {S,H}=geopar; \
{Em}=matpar; {A10,A20,A30}=fabpar; \n        A0=(A10+A20)/2;  {uX,uY}=state; \
\n        q=IncLoadOfMisesTruss[{S,H,Em,A0,0},\[Lambda],{uX,uY},numer];\n     \
   If [!numer,q=Simplify[q]]; Return[q]];\n    If [pnam2==\"CG\", \
\[Mu]=state[[1]];  bipath=False;\n        If [pnam3==\"CGB\", bipath=True]; \n\
        q=IncLoadOfCircleGame[{bipath},\[Lambda],{\[Mu]},numer];\n        If \
[!numer,q=Simplify[q]]; Return[q]];\n    Print[\" Problem \",problem,\" not \
implemented\"];\n    Return[Null]];\n",
  StyleBox["    \n\
(*ClearAll[S,H,Em,A0,uX,uY,\[Lambda],\[Mu],numer,ptags,geopar,\n\
matpar,fabpar,solpar,",
    FontColor->RGBColor[0, 0, 1]],
  StyleBox["reffor",
    FontColor->RGBColor[0, 0, 1]],
  StyleBox[",conpar,state,options];\nptags={\"MT\",\"FE\",\"LC\"}; \
geopar={S,H}; matpar={Em}; fabpar={A0,A0,0};  \nsolpar={}; \
conpar={\[Lambda]}; state={uX,uY}; numer=False; options={numer};\n\
q=IncLoad[ptags,geopar,matpar,fabpar,solpar,\n          ",
    FontColor->RGBColor[0, 0, 1]],
  StyleBox["reffor",
    FontColor->RGBColor[0, 0, 1]],
  StyleBox[",conpar,state,options]; Print[\"q=\",q//MatrixForm];\n\
ptags={\"CGB\",\"FE\",\"LC\"}; state={\[Mu]}; \n\
q=IncLoad[ptags,geopar,matpar,fabpar,solpar,\n          ",
    FontColor->RGBColor[0, 0, 1]],
  StyleBox["reffor",
    FontColor->RGBColor[0, 0, 1]],
  StyleBox[",conpar,state,options]; Print[\"q=\",q//MatrixForm];*)",
    FontColor->RGBColor[0, 0, 1]]
}], "Input",
  CellFrame->True,
  CellMargins->{{14, 37}, {Inherited, Inherited}},
  InitializationCell->True,
  CellLabelMargins->{{8, Inherited}, {Inherited, Inherited}},
  ImageRegion->{{-0, 1}, {0, 1}},
  Background->RGBColor[0, 1, 0]],

Cell[TextData[{
  StyleBox["Solution advancing modules",
    FontWeight->"Bold"],
  StyleBox["  based on purely-incremental continuation",
    FontVariations->{"CompatibilityType"->0}]
}], "Text",
  CellFrame->True,
  CellMargins->{{14, 37}, {Inherited, Inherited}},
  CellLabelMargins->{{8, Inherited}, {Inherited, Inherited}},
  ImageRegion->{{-0, 1}, {0, 1}},
  Background->RGBColor[0, 1, 1]],

Cell[TextData[{
  StyleBox["IncVel",
    FontWeight->"Bold"],
  StyleBox[":  ",
    FontWeight->"Bold"],
  StyleBox["computes incremental velocity vector",
    FontVariations->{"CompatibilityType"->0}],
  " (RHS of incremental ODE).\nPerforms singularity and ill-condition tests."
}], "Text",
  CellFrame->True,
  CellMargins->{{14, 37}, {Inherited, Inherited}},
  CellLabelMargins->{{8, Inherited}, {Inherited, Inherited}},
  ImageRegion->{{-0, 1}, {0, 1}},
  Background->RGBColor[1, 1, 0]],

Cell[TextData[{
  "IncVel[ptags_,geopar_,matpar_,fabpar_,solpar_,reffor_,\n  \
conpar_,state_,options_]:= Module[{K,Kdet,Knorm,\n  \
Knill,Kcinv,numer,ev,evmin,evmax,Null3=Table[Null,{3}],\n  \
problem,int,ics,epsill=10.^(-6),epsing=10.^(-12),\n  status=\" \",q,v}, \
{problem,int,ics}=ptags;\n",
  StyleBox[" (*Print[\"Enter IncVel, problem=\",problem,\"  \
method=\",method];*)\n (* Print[\"geopar=\",geopar,\" matpar=\",matpar,\" \
fabpar=\",fabpar];*)\n (*Print[\"solpar=\",solpar];*)\n (* \
Print[\"conpar=\",conpar,\" state=\",state];*)",
    FontColor->RGBColor[1, 0, 0]],
  "\n  K=TanStiff[ptags,geopar,matpar,fabpar,solpar,\n             \
reffor,conpar,state,options];",
  StyleBox[" \n  (*Print[\"K=\",K//MatrixForm];*)",
    FontColor->RGBColor[1, 0, 0]],
  "\n  If [K==Null, status=\"Error in K eval\"; Return[{Null3,status}]];\n  \
Kdet=Det[K]; numer=options[[1]];\n  If [!numer, Kdet=Simplify[Kdet];",
  StyleBox[" ",
    FontColor->RGBColor[1, 0, 0]],
  "\n      If [Kdet==0, status=\"Singular K\"; Return[{Null3,status}]]];\n  \
If [numer, Knorm=MatrixFrobNorm[K]; Knill=epsill*Knorm;\n      If [Knorm==0, \
status=\"Singular K\"; Return[{Null3,status}]];\n",
  StyleBox["      (*Print[\"Kdet,Knorm=\",{Kdet,Knorm},\" Knill=\",Knill];*)\n\
      (*Print[\"check Abs[Kdet]<=Knill= \",Abs[Kdet]<=Knill]; *)",
    FontColor->RGBColor[1, 0, 0]],
  "\n      If [Abs[Kdet]<=Knill, ev=Abs[Eigenvalues[N[K]]]; ",
  StyleBox["\n",
    FontColor->RGBColor[1, 0, 0]],
  "          evmax=Max[Max[ev],epsing]; evmin=Min[ev];  \n        ",
  StyleBox["  (*Print[\"ev=\",ev,\"  evmin,max=\",{evmin,evmax}];*)",
    FontColor->RGBColor[1, 0, 0]],
  "\n          Kcinv=evmin/evmax; status=\"Ill conditioned K\";  \n        ",
  StyleBox["  (*Print[\"Kcinv=\",Kcinv];*)",
    FontColor->RGBColor[1, 0, 0]],
  "\n          If [Kcinv<=epsing, status=\"Singular K\"; \n              \
Return[{Null3,status}]]];\n      ];       \n  \
q=IncLoad[ptags,geopar,matpar,fabpar,solpar,\n            \
reffor,conpar,state,options];",
  StyleBox[" \n  (*Print[\"q=\",q//MatrixForm];*)",
    FontColor->RGBColor[1, 0, 0]],
  "\n  v=LinearSolve[K,q];   If [!numer,v=Simplify[v]]; \n",
  StyleBox["  (*Print[\"after LinearSolve, K=\",K//MatrixForm]; *)\n  \
(*Print[\"v=\",v//MatrixForm]; Print[\"status=\",status];*)",
    FontColor->RGBColor[1, 0, 0]],
  "\n  Return[{{v,q,Kdet},status}]];\n",
  StyleBox["\n(*ClearAll[S,H,Em,A0,uX,uY,\[Lambda],\[Mu],numer,ptags,geopar,\n\
matpar,fabpar,solpar,",
    FontColor->RGBColor[0, 0, 1]],
  StyleBox["reffor",
    FontColor->RGBColor[0, 0, 1]],
  StyleBox[",conpar,state,options]; \nptags={\"MT\",\"FE\",\"LC\"}; \
numer=True; \nS=2; H=1; Em=10; A0=1; uX=0; uY=(Sqrt[3]-3)/3+.000000001;\n\
geopar={S,H}; matpar={Em}; fabpar={A0,A0,0};  \nsolpar={}; \
conpar={\[Lambda]}; state={uX,uY}; options={numer};\nIf [numer, \
{geopar,matpar,fabpar,state}=N[{geopar,matpar,fabpar,state}]]; \n\
{{v,q,Kdet},status}=IncVel[ptags,geopar,matpar,fabpar,solpar,\n               \
   ",
    FontColor->RGBColor[0, 0, 1]],
  StyleBox["reffor",
    FontColor->RGBColor[0, 0, 1]],
  StyleBox[",conpar,state,options]; \nPrint[\"v=\",v//MatrixForm,\" \
q=\",q//MatrixForm,\" Kdet=\",Kdet]; \nIf [status!=\" \",Print[\"status: \
\",status]];\nptags={\"CGB\",\"FE\",\"LC\"}; state={\[Mu]}; \[Mu]=2/3; \
\[Lambda]=1; conpar={\[Lambda]}; \nnumer=False; options={numer};\n\
{{v,q,Kdet},status}=IncVel[ptags,geopar,matpar,fabpar,solpar,\n               \
   ",
    FontColor->RGBColor[0, 0, 1]],
  StyleBox["reffor",
    FontColor->RGBColor[0, 0, 1]],
  StyleBox[",conpar,state,options]; \nPrint[\"v=\",v//MatrixForm,\" \
q=\",q//MatrixForm,\" Kdet=\",Kdet]; \nIf [status!=\" \",Print[\"status: \
\",status]];*)",
    FontColor->RGBColor[0, 0, 1]]
}], "Input",
  CellFrame->True,
  CellMargins->{{14, 37}, {Inherited, Inherited}},
  InitializationCell->True,
  CellLabelMargins->{{8, Inherited}, {Inherited, Inherited}},
  ImageRegion->{{-0, 1}, {0, 1}},
  Background->RGBColor[0, 1, 0]],

Cell[TextData[{
  StyleBox["IncStep1: ",
    FontWeight->"Bold",
    FontColor->GrayLevel[0]],
  StyleBox["step increment module for a one-stage, explicit RK integrator.  \
Since Forward \nEuler (FE) is the only instance of this class, a call to ",
    FontColor->GrayLevel[0]],
  StyleBox["ButcherTable",
    FontWeight->"Bold",
    FontColor->GrayLevel[0]],
  StyleBox[" is nor needed",
    FontColor->GrayLevel[0]]
}], "Text",
  CellFrame->True,
  CellMargins->{{14, 37}, {Inherited, Inherited}},
  CellLabelMargins->{{8, Inherited}, {Inherited, Inherited}},
  ImageRegion->{{-0, 1}, {0, 1}},
  Background->RGBColor[1, 1, 0]],

Cell[TextData[{
  StyleBox["IncStep1[",
    FontColor->GrayLevel[0]],
  "ptags_",
  StyleBox[",geopar_,matpar_,fabpar_,solpar_,",
    FontColor->GrayLevel[0]],
  "reffor_,\n   conpar_,state_,options_,soldim_",
  StyleBox["]:=Module[\n  {nmax,\[Lambda]max,umax,elln,ellmin,ellmax,acctol,\n\
   problem,",
    FontColor->GrayLevel[0]],
  "int,ics,cscale,nextsol",
  StyleBox[",ell,a=1,\[Kappa]=1,n,m,\n   un,v1,vw,f1,fw,q1,qw,Kdet1,",
    FontColor->GrayLevel[0]],
  "\[Lambda]n",
  StyleBox[",",
    FontColor->GrayLevel[0]],
  "\[CapitalDelta]\[Lambda]n,unext,\[Lambda]next",
  StyleBox["},\n  ",
    FontColor->GrayLevel[0]],
  "{problem,int,ics}=ptags;",
  StyleBox[" \n ",
    FontColor->GrayLevel[0]],
  StyleBox["  (*Print[\"Enter IncStep1\"];*)",
    FontColor->RGBColor[1, 0, 0]],
  StyleBox["\n  {nmax,\[Lambda]max,umax,ell,ellmin,ellmax,acctol,",
    FontColor->GrayLevel[0]],
  "cscale}",
  StyleBox["=solpar;\n  ",
    FontColor->GrayLevel[0]],
  StyleBox[" (*Print[\"solpar=\",solpar];*)",
    FontColor->RGBColor[1, 0, 0]],
  StyleBox["\n",
    FontColor->GrayLevel[0]],
  StyleBox["   \
(*Print[\"nmax,\[Lambda]max,umax,ell,ellmin,ellmax,acctol,cscale=\",\n        \
 {nmax,\[Lambda]max,umax,ell,ellmin,ellmax,acctol,cscale}];*)",
    FontColor->RGBColor[1, 0, 0]],
  StyleBox["\n   {n,m}=soldim; ",
    FontColor->GrayLevel[0]],
  "\[Lambda]n",
  StyleBox["=conpar[[1]]; un=state; \n",
    FontColor->GrayLevel[0]],
  StyleBox["   (*Print[\"n=\",n,\" \[Lambda]n=\",\[Lambda]n];*)",
    FontColor->RGBColor[1, 0, 0]],
  StyleBox["\n   nextsol=Table[Null,{m}]; ",
    FontColor->GrayLevel[0]],
  "elln=ell; ",
  StyleBox["\n  {{v1,q1,Kdet1},status}=",
    FontColor->GrayLevel[0]],
  "IncVel[ptags,",
  StyleBox["geopar,matpar,fabpar,\n                 solpar,",
    FontColor->GrayLevel[0]],
  "reffor,conpar,state,options",
  StyleBox["];\n",
    FontColor->GrayLevel[0]],
  StyleBox["   (*Print[\"v1=\",v1//MatrixForm,\", q1=\",q1//MatrixForm,\n     \
  \", Kdet1=\",Kdet1];*)",
    FontColor->RGBColor[1, 0, 0]],
  StyleBox["\n   If [status!=\" \",Return[{nextsol,status}]];",
    FontColor->GrayLevel[0]],
  "\n   f1=ICSFactor[ics,cscale,v1,q1]; \[CapitalDelta]\[Lambda]1=elln/f1; \n \
",
  StyleBox["  (*Print[\"f1=\",f1,\", \[CapitalDelta]\[Lambda]1=\",\
\[CapitalDelta]\[Lambda]1];*)",
    FontColor->RGBColor[1, 0, 0]],
  "\n   vw=v1; qw=q1; fw=f1;\n   \[CapitalDelta]\[Lambda]n=elln/fw;  \
\[Lambda]next=\[Lambda]n+\[CapitalDelta]\[Lambda]n; unext=un+vw*\
\[CapitalDelta]\[Lambda]n; n++;\n ",
  StyleBox["  \
(*Print[\"\[CapitalDelta]\[Lambda]n=\",\[CapitalDelta]\[Lambda]n,\" \
\[Lambda]next=\",\[Lambda]next,\", unext=\",unext];*)",
    FontColor->RGBColor[1, 0, 0]],
  "\n   rnext=TotalResidual[ptags,geopar,matpar,fabpar,\n         \
solpar,reffor,{\[Lambda]next",
  StyleBox["},unext,",
    FontColor->GrayLevel[0]],
  "options];",
  StyleBox[" \n   (*Print[\"rnext=\",rnext];*)",
    FontColor->RGBColor[1, 0, 0]],
  "\n   nextsol",
  StyleBox["={n,\[Lambda]next,unext,vw,rnext,qw,Kdet1,ell,a,\[Kappa],\" \"};\n\
",
    FontColor->GrayLevel[0]],
  StyleBox["   (*Print[\"IncStep1 nextsol=\",nextsol];*)",
    FontColor->RGBColor[1, 0, 0]],
  StyleBox["\n   Return[{",
    FontColor->GrayLevel[0]],
  "nextsol",
  StyleBox[",\" \"}]];",
    FontColor->GrayLevel[0]],
  "\n",
  StyleBox["   \n",
    FontColor->RGBColor[0, 0, 1]],
  StyleBox["(*ClearAll[S,H,Em,A0,uX,uY,\[Lambda],\[Mu],numer,ptags,geopar,\
matpar,\nfabpar,solpar,reffor,conpar,state,options,incr,cscale];\nptags={\"MT\
\",\"EMR\",\"LC\"}; numer=True;  incr=0;\nS=2; H=1; Em=10; A0=1; \[Lambda]=0; \
 uX=0; uY=(Sqrt[3]-3)/3; uY=0;\ngeopar={S,H}; matpar={Em}; fabpar={A0,A0,0};\n\
nmax=1; \[Lambda]max=1; umax=1; ell=1; a=0; \[Kappa]=0; \[Kappa]0=1; \n\
ellmin=.1; ellmax=10.; acctol=.001; cscale={}; rem=\" \";\nsolpar={nmax,\
\[Lambda]max,umax,ell,ellmin,ellmax,acctol,cscale}; \nconpar={\[Lambda]}; \
state={uX,uY}; options={numer}; soldim={0,11};\nIf [numer, \
{geopar,matpar,fabpar,state}=N[{geopar,matpar,fabpar,state}]]; \n\
{nextsol,status}=IncStep1[ptags,geopar,matpar,fabpar,solpar,\n                \
 reffor,conpar,state,options,soldim];\nPrint[\"nextsol=\",nextsol]; \n\
If[status!=\" \",Print[\"status: \",status]];*)",
    FontColor->RGBColor[0, 0, 1]]
}], "Input",
  CellFrame->True,
  CellMargins->{{14, 37}, {Inherited, Inherited}},
  InitializationCell->True,
  CellLabelMargins->{{8, Inherited}, {Inherited, Inherited}},
  ImageRegion->{{-0, 1}, {0, 1}},
  Background->RGBColor[0, 1, 0]],

Cell[TextData[{
  StyleBox["IncStep2: ",
    FontWeight->"Bold",
    FontColor->GrayLevel[0]],
  StyleBox["step increment module for a two-stage, explicit, second-order RK \
integrator. \nThis is a one-parameter family that includes Explicit Midpoint \
Rule (EMR) and Explicit\nTrapezoidal Rule (ETR) as key instances. ",
    FontColor->GrayLevel[0]],
  StyleBox["ButcherTable",
    FontWeight->"Bold",
    FontColor->GrayLevel[0],
    FontVariations->{"CompatibilityType"->0}],
  StyleBox[" supplies method coefficients",
    FontColor->GrayLevel[0]]
}], "Text",
  CellFrame->True,
  CellMargins->{{14, 37}, {Inherited, Inherited}},
  CellLabelMargins->{{8, Inherited}, {Inherited, Inherited}},
  ImageRegion->{{-0, 1}, {0, 1}},
  Background->RGBColor[1, 1, 0]],

Cell[TextData[{
  StyleBox["IncStep2[",
    FontColor->GrayLevel[0]],
  "ptags_",
  StyleBox[",geopar_,matpar_,fabpar_,solpar_,",
    FontColor->GrayLevel[0]],
  "reffor_,\n   conpar_,state_,options_,soldim_",
  StyleBox["]:=Module[\n  {nmax,\[Lambda]max,umax,elln,ellmin,ellmax,acctol,",
    
    FontColor->GrayLevel[0]],
  "cscale,",
  StyleBox["\n   problem,",
    FontColor->GrayLevel[0]],
  "int,ics,nextsol",
  StyleBox[",ell,a=1,\[Kappa]=1,A,c,w,c1,c2,\n   \
w1,w2,A21,n,m,un,u1,u2,v1,v2,vw,q1,q2,qw,",
    FontColor->GrayLevel[0]],
  "\[Lambda]n",
  StyleBox[",",
    FontColor->GrayLevel[0]],
  "\[Lambda]1,\[Lambda]2,\n   \[CapitalDelta]\[Lambda]n,\[CapitalDelta]\
\[Lambda]1,\[CapitalDelta]\[Lambda]2,f1,f2,fw,",
  StyleBox["Kdet1,Kdet2,",
    FontColor->GrayLevel[0]],
  "unext,\[Lambda]next",
  StyleBox["}, \n  ",
    FontColor->GrayLevel[0]],
  "{problem,int,ics}=ptags;",
  StyleBox[" \n  {nmax,\[Lambda]max,umax,ell,ellmin,ellmax,acctol,",
    FontColor->GrayLevel[0]],
  "cscale}",
  StyleBox["=solpar;\n  {n,m}=soldim; ",
    FontColor->GrayLevel[0]],
  "\[Lambda]n",
  StyleBox["=conpar[[1]]; un=state;\n   nextsol=Table[Null,{m}]; ",
    FontColor->GrayLevel[0]],
  "elln=ell;\n  {A,c,w}=ButcherTable[int]; \n  {c1,c2}=c; {w1,w2}=w; \
A21=A[[2,1]]; ",
  StyleBox["\n  {{v1,q1,Kdet1},status}=",
    FontColor->GrayLevel[0]],
  "IncVel[ptags,",
  StyleBox["\n    geopar,matpar,fabpar,solpar,",
    FontColor->GrayLevel[0]],
  "reffor",
  StyleBox[",{\[Lambda]n},un,options];\n   If [status!=\" \",",
    FontColor->GrayLevel[0]],
  "Return[{nextsol,status}]",
  StyleBox["];",
    FontColor->GrayLevel[0]],
  "\n   f1=ICSFactor[ics,cscale,v1,q1]; \[CapitalDelta]\[Lambda]1=elln/f1; \n \
  \[Lambda]2=\[Lambda]n+c2*\[CapitalDelta]\[Lambda]1; u2=un+A21*v1*\
\[CapitalDelta]\[Lambda]1;\n  ",
  StyleBox["{{v2,q2,Kdet2},status}=",
    FontColor->GrayLevel[0]],
  "IncVel[ptags,",
  StyleBox["\n    geopar,matpar,fabpar,solpar,",
    FontColor->GrayLevel[0]],
  "reffor",
  StyleBox[",{\[Lambda]2},u2,options];\n   If [status!=\" \",",
    FontColor->GrayLevel[0]],
  "Return[{nextsol,status}]",
  StyleBox["];",
    FontColor->GrayLevel[0]],
  "\n   vw=w1*v1+w2*v2;  qw=w1*q1+w2*q2;\n   fw=ICSFactor[ics,cscale,vw,qw];\n\
   \[CapitalDelta]\[Lambda]n=elln/fw;  \[Lambda]next=\[Lambda]n+\
\[CapitalDelta]\[Lambda]n; unext=un+vw*\[CapitalDelta]\[Lambda]n; n++;\n   \
rnext=TotalResidual[ptags,geopar,matpar,fabpar,\n         solpar,reffor,{\
\[Lambda]next",
  StyleBox["},unext,",
    FontColor->GrayLevel[0]],
  "options];",
  StyleBox[" ",
    FontColor->RGBColor[1, 0, 0]],
  "\n   nextsol",
  StyleBox["={n,\[Lambda]next,unext,vw,rnext,qw,Kdet2,ell,a,\[Kappa],\" \"};\n\
   Return[{",
    FontColor->GrayLevel[0]],
  "nextsol",
  StyleBox[",\" \"}]];",
    FontColor->GrayLevel[0]],
  " "
}], "Input",
  CellFrame->True,
  CellMargins->{{14, 37}, {Inherited, Inherited}},
  InitializationCell->True,
  CellLabelMargins->{{8, Inherited}, {Inherited, Inherited}},
  ImageRegion->{{-0, 1}, {0, 1}},
  Background->RGBColor[0, 1, 0]],

Cell[TextData[{
  StyleBox["IncStep3: ",
    FontWeight->"Bold",
    FontColor->GrayLevel[0]],
  StyleBox["step increment module for a three-stage, explicit, third-order RK \
integrator. \nThis is a multiparameter family that includes Ralston's version \
(RK3R) and Heun's\nversion (RK3H) as key instances.  ",
    FontColor->GrayLevel[0]],
  StyleBox["ButcherTable",
    FontWeight->"Bold",
    FontColor->GrayLevel[0],
    FontVariations->{"CompatibilityType"->0}],
  StyleBox[" supplies method coefficients",
    FontColor->GrayLevel[0]]
}], "Text",
  CellFrame->True,
  CellMargins->{{14, 37}, {Inherited, Inherited}},
  CellLabelMargins->{{8, Inherited}, {Inherited, Inherited}},
  ImageRegion->{{-0, 1}, {0, 1}},
  Background->RGBColor[1, 1, 0]],

Cell[TextData[{
  StyleBox["IncStep3[",
    FontColor->GrayLevel[0]],
  "ptags_",
  StyleBox[",geopar_,matpar_,fabpar_,solpar_,",
    FontColor->GrayLevel[0]],
  "reffor_,\n   conpar_,state_,options_,soldim_",
  StyleBox["]:=Module[\n  {nmax,\[Lambda]max,umax,elln,ellmin,ellmax,acctol,",
    
    FontColor->GrayLevel[0]],
  "cscale,",
  StyleBox["\n   problem,",
    FontColor->GrayLevel[0]],
  "int,ics,nextsol",
  StyleBox[",ell,a=1,\[Kappa]=1,A,c,w,c1,c2,c3,\n   \
w1,w2,w3,A21,A31,A32,n,m,un,u1,u2,u3,v1,v2,v3,vw,\n   q1,q2,q3,qw,",
    FontColor->GrayLevel[0]],
  "\[Lambda]n",
  StyleBox[",",
    FontColor->GrayLevel[0]],
  "\[Lambda]1,\[Lambda]2,\[Lambda]3,\[CapitalDelta]\[Lambda]n,\[CapitalDelta]\
\[Lambda]1,\[CapitalDelta]\[Lambda]2,\[CapitalDelta]\[Lambda]3,f1,f2,f3,fw,",
  StyleBox["\n   Kdet1,Kdet2,Kdet3,",
    FontColor->GrayLevel[0]],
  "unext,\[Lambda]next",
  StyleBox["}, \n  ",
    FontColor->GrayLevel[0]],
  "{problem,int,ics}=ptags;",
  StyleBox[" \n  {nmax,\[Lambda]max,umax,ell,ellmin,ellmax,acctol,",
    FontColor->GrayLevel[0]],
  "cscale}",
  StyleBox["=solpar;\n  {n,m}=soldim;; ",
    FontColor->GrayLevel[0]],
  "\[Lambda]n",
  StyleBox["=conpar[[1]]; un=state; \n   nextsol=Table[Null,{m}]; ",
    FontColor->GrayLevel[0]],
  "elln=ell;\n  {A,c,w}=ButcherTable[int]; {c1,c2,c3}=c; \n  {w1,w2,w3}=w; \
A21=A[[2,1]]; A31=A[[3,1]]; A32=A[[3,2]];",
  StyleBox["\n  {{v1,q1,Kdet1},status}=",
    FontColor->GrayLevel[0]],
  "IncVel[ptags,",
  StyleBox["\n        geopar,matpar,",
    FontColor->GrayLevel[0]],
  "fabpar,solpar,reffor,{\[Lambda]",
  StyleBox["n},un,options];\n   If [status!=\" \",",
    FontColor->GrayLevel[0]],
  "Return[{nextsol,status}]",
  StyleBox["];",
    FontColor->GrayLevel[0]],
  "\n   f1=ICSFactor[ics,cscale,v1,q1]; \[CapitalDelta]\[Lambda]1=elln/f1; \n \
  \[Lambda]2=\[Lambda]n+c2*\[CapitalDelta]\[Lambda]1; u2=un+A21*v1*\
\[CapitalDelta]\[Lambda]1;\n  ",
  StyleBox["{{v2,q2,Kdet2},status}=",
    FontColor->GrayLevel[0]],
  "IncVel[ptags,",
  StyleBox["\n        geopar,matpar,",
    FontColor->GrayLevel[0]],
  "fabpar,solpar,reffor,{\[Lambda]",
  StyleBox["2},u2,options];\n   If [status!=\" \",",
    FontColor->GrayLevel[0]],
  "Return[{nextsol,status}]",
  StyleBox["];\n   ",
    FontColor->GrayLevel[0]],
  "f2=ICSFactor[ics,cscale,v2,q2]; \[CapitalDelta]\[Lambda]2=elln/f2; \n   \
\[Lambda]3=\[Lambda]n+c3*\[CapitalDelta]\[Lambda]2; u3=un+(A31*v1+A32*v2)*\
\[CapitalDelta]\[Lambda]2;\n  ",
  StyleBox["{{v3,q3,Kdet3},status}=",
    FontColor->GrayLevel[0]],
  "IncVel[ptags,",
  StyleBox["\n        geopar,matpar,",
    FontColor->GrayLevel[0]],
  "fabpar,solpar,reffor,{\[Lambda]",
  StyleBox["3},u3,options];\n   If [status!=\" \",",
    FontColor->GrayLevel[0]],
  "Return[{nextsol,status}]",
  StyleBox["];",
    FontColor->GrayLevel[0]],
  "\n   vw=w1*v1+w2*v2+w3*v3;  qw=w1*q1+w2*q2+w3*q3;\n   \
fw=ICSFactor[ics,cscale,vw,qw];\n   \[CapitalDelta]\[Lambda]n=elln/fw;  \
\[Lambda]next=\[Lambda]n+\[CapitalDelta]\[Lambda]n; unext=un+vw*\
\[CapitalDelta]\[Lambda]n; n++;\n   \
rnext=TotalResidual[ptags,geopar,matpar,fabpar,\n         solpar,reffor,{\
\[Lambda]next",
  StyleBox["},unext,",
    FontColor->GrayLevel[0]],
  "options];",
  StyleBox[" ",
    FontColor->RGBColor[1, 0, 0]],
  "\n   nextsol",
  StyleBox["={n,\[Lambda]next,unext,vw,rnext,qw,Kdet3,ell,a,\[Kappa],\" \"};\n\
   Return[{",
    FontColor->GrayLevel[0]],
  "nextsol",
  StyleBox[",\" \"}]];",
    FontColor->GrayLevel[0]],
  " "
}], "Input",
  CellFrame->True,
  CellMargins->{{14, 37}, {Inherited, Inherited}},
  InitializationCell->True,
  CellLabelMargins->{{8, Inherited}, {Inherited, Inherited}},
  ImageRegion->{{-0, 1}, {0, 1}},
  Background->RGBColor[0, 1, 0]],

Cell[TextData[{
  StyleBox["IncStep4: ",
    FontWeight->"Bold",
    FontColor->GrayLevel[0]],
  StyleBox["step increment module for a four-stage, explicit, fourth-order RK \
integrator. \nThis is a multiparameter family that includes the classic \
version (RK4C) and Kutta's\nversion (RK4K) as key instances. ",
    FontColor->GrayLevel[0]],
  StyleBox["ButcherTable",
    FontWeight->"Bold",
    FontColor->GrayLevel[0],
    FontVariations->{"CompatibilityType"->0}],
  StyleBox[" supplies method coefficients",
    FontColor->GrayLevel[0]]
}], "Text",
  CellFrame->True,
  CellMargins->{{14, 37}, {Inherited, Inherited}},
  CellLabelMargins->{{8, Inherited}, {Inherited, Inherited}},
  ImageRegion->{{-0, 1}, {0, 1}},
  Background->RGBColor[1, 1, 0]],

Cell[TextData[{
  StyleBox["IncStep4[",
    FontColor->GrayLevel[0]],
  "ptags_",
  StyleBox[",geopar_,matpar_,fabpar_,solpar_,",
    FontColor->GrayLevel[0]],
  "reffor_,\n   conpar_,state_,options_,soldim_",
  StyleBox["]:=Module[\n  {nmax,\[Lambda]max,umax,elln,ellmin,ellmax,acctol,",
    
    FontColor->GrayLevel[0]],
  "cscale,\n   ",
  StyleBox["problem,",
    FontColor->GrayLevel[0]],
  "int,ics,nextsol",
  StyleBox[",ell,a=1,\[Kappa]=1,A,c,w,c1,c2,c3,c4,\n   \
w1,w2,w3,w4,A21,A31,A32,A41,A42,A43,\n   \
n,m,un,u1,u2,u3,u4,v1,v2,v3,v4,vw,q1,q2,q3,q4,qw,\n   ",
    FontColor->GrayLevel[0]],
  "\[Lambda]n",
  StyleBox[",",
    FontColor->GrayLevel[0]],
  "\[Lambda]1,\[Lambda]2,\[Lambda]3,\[Lambda]4,\[CapitalDelta]\[Lambda]n,\
\[CapitalDelta]\[Lambda]1,\[CapitalDelta]\[Lambda]2,\[CapitalDelta]\[Lambda]3,\
\[CapitalDelta]\[Lambda]4,f1,f2,f3,f4,fw,",
  StyleBox["\n   Kdet1,Kdet2,Kdet3,Kdet4,",
    FontColor->GrayLevel[0]],
  "unext,\[Lambda]next",
  StyleBox["}, \n  ",
    FontColor->GrayLevel[0]],
  "{problem,int,ics}=ptags;",
  StyleBox[" \n  {nmax,\[Lambda]max,umax,ell,ellmin,ellmax,acctol,",
    FontColor->GrayLevel[0]],
  "cscale}",
  StyleBox["=solpar;\n  {n,m}=soldim; ",
    FontColor->GrayLevel[0]],
  "\[Lambda]n",
  StyleBox["=conpar[[1]]; un=state; \n   nextsol=Table[Null,{m}]; ",
    FontColor->GrayLevel[0]],
  "elln=ell; \n  {A,c,w}=ButcherTable[int]; {c1,c2,c3,c4}=c; \n  \
{w1,w2,w3,w4}=w; A21=A[[2,1]]; A31=A[[3,1]]; \n   A32=A[[3,2]]; A41=A[[4,1]]; \
A42=A[[4,2]]; A43=A[[4,3]]; ",
  StyleBox["\n  {{v1,q1,Kdet1},status}=",
    FontColor->GrayLevel[0]],
  "IncVel[ptags,",
  StyleBox["\n    geopar,matpar,",
    FontColor->GrayLevel[0]],
  "fabpar,solpar,reffor,{\[Lambda]",
  StyleBox["n},un,options];\n   If [status!=\" \",",
    FontColor->GrayLevel[0]],
  "Return[{nextsol,status}]",
  StyleBox["];",
    FontColor->GrayLevel[0]],
  "\n   f1=ICSFactor[ics,cscale,v1,q1]; \[CapitalDelta]\[Lambda]1=elln/f1; \n \
  \[Lambda]2=\[Lambda]n+c2*\[CapitalDelta]\[Lambda]1; u2=un+A21*v1*\
\[CapitalDelta]\[Lambda]1;\n  ",
  StyleBox["{{v2,q2,Kdet2},status}=",
    FontColor->GrayLevel[0]],
  "IncVel[ptags,",
  StyleBox["\n     geopar,matpar,",
    FontColor->GrayLevel[0]],
  "fabpar,solpar,reffor,{\[Lambda]",
  StyleBox["2},u2,options];\n   If [status!=\" \",",
    FontColor->GrayLevel[0]],
  "Return[{nextsol,status}]",
  StyleBox["];\n   ",
    FontColor->GrayLevel[0]],
  "f2=ICSFactor[ics,cscale,v2,q2]; \[CapitalDelta]\[Lambda]2=elln/f2; \n   \
\[Lambda]3=\[Lambda]n+c3*\[CapitalDelta]\[Lambda]2; u3=un+(A31*v1+A32*v2)*\
\[CapitalDelta]\[Lambda]2;\n  ",
  StyleBox["{{v3,q3,Kdet3},status}=",
    FontColor->GrayLevel[0]],
  "IncVel[ptags,",
  StyleBox["\n     geopar,matpar,",
    FontColor->GrayLevel[0]],
  "fabpar,solpar,reffor,{\[Lambda]",
  StyleBox["3},u3,options];\n   If [status!=\" \",",
    FontColor->GrayLevel[0]],
  "Return[{nextsol,status}]",
  StyleBox["];\n   ",
    FontColor->GrayLevel[0]],
  "f3=ICSFactor[ics,cscale,v3,q3]; \[CapitalDelta]\[Lambda]3=elln/f3; \n   \
\[Lambda]4=\[Lambda]n+c4*\[CapitalDelta]\[Lambda]3; \
u4=un+(A41*v1+A42*v2+A43*v3)*\[CapitalDelta]\[Lambda]3;\n  ",
  StyleBox["{{v4,q4,Kdet4},status}=",
    FontColor->GrayLevel[0]],
  "IncVel[ptags,",
  StyleBox["\n    geopar,matpar,",
    FontColor->GrayLevel[0]],
  "fabpar,solpar,reffor,{\[Lambda]",
  StyleBox["4},u4,options];\n   If [status!=\" \",",
    FontColor->GrayLevel[0]],
  "Return[{nextsol,status}]",
  StyleBox["];",
    FontColor->GrayLevel[0]],
  "\n   vw=w1*v1+w2*v2+w3*v3+w4*v4; qw=w1*q1+w2*q2+w3*q3+w4*q4;\n   \
fw=ICSFactor[ics,cscale,vw,qw];\n   \[CapitalDelta]\[Lambda]n=elln/fw;  \
\[Lambda]next=\[Lambda]n+\[CapitalDelta]\[Lambda]n; unext=un+vw*\
\[CapitalDelta]\[Lambda]n; n++;\n   \
rnext=TotalResidual[ptags,geopar,matpar,fabpar,\n         solpar,reffor,{\
\[Lambda]next",
  StyleBox["},unext,",
    FontColor->GrayLevel[0]],
  "options];",
  StyleBox[" ",
    FontColor->RGBColor[1, 0, 0]],
  "\n   nextsol",
  StyleBox["={n,\[Lambda]next,unext,vw,rnext,qw,Kdet4,ell,a,\[Kappa],\" \"};\n\
   Return[{",
    FontColor->GrayLevel[0]],
  "nextsol",
  StyleBox[",\" \"}]];",
    FontColor->GrayLevel[0]],
  " "
}], "Input",
  CellFrame->True,
  CellMargins->{{14, 37}, {Inherited, Inherited}},
  InitializationCell->True,
  CellLabelMargins->{{8, Inherited}, {Inherited, Inherited}},
  ImageRegion->{{-0, 1}, {0, 1}},
  Background->RGBColor[0, 1, 0]],

Cell[TextData[{
  StyleBox["IncStep",
    FontWeight->"Bold"],
  ": step advancing cover module.  Branches to one of ",
  StyleBox["IncStep1",
    FontWeight->"Bold"],
  " through ",
  StyleBox["IncStep4 ",
    FontWeight->"Bold"],
  "as per \nnumber of RK evaluations (aka stages) over step.  This # is \
checked via ",
  StyleBox["ButcherTable",
    FontWeight->"Bold"],
  ". "
}], "Text",
  CellFrame->True,
  CellMargins->{{14, 37}, {Inherited, Inherited}},
  CellLabelMargins->{{8, Inherited}, {Inherited, Inherited}},
  ImageRegion->{{-0, 1}, {0, 1}},
  Background->RGBColor[1, 1, 0]],

Cell[TextData[{
  "IncStep",
  StyleBox["[",
    FontColor->GrayLevel[0]],
  "ptags_",
  StyleBox[",geopar_,matpar_,fabpar_,solpar_,",
    FontColor->GrayLevel[0]],
  "reffor_,\n        conpar_,state_,options_,soldim_",
  StyleBox["]:=Module[\n  {problem,int,ics,A,c,w,neval,nextsol,status},\n  ",
    FontColor->GrayLevel[0]],
  StyleBox["(*Print[\"Enter IncStep\"];*)",
    FontColor->RGBColor[1, 0, 0]],
  StyleBox["\n  ",
    FontColor->GrayLevel[0]],
  "{problem,int,ics}=ptags;",
  StyleBox[" ",
    FontColor->GrayLevel[0]],
  "{A,c,w}=ButcherTable[int]; \n   neval=Length[c];\n   If [neval==1,\n     \
{nextsol,status}=IncStep1[ptags,geopar,matpar,fabpar,\n      \
solpar,reffor,conpar,state,options,",
  StyleBox["soldim]",
    FontColor->GrayLevel[0]],
  "];\n   If [neval==2,\n     \
{nextsol,status}=IncStep2[ptags,geopar,matpar,fabpar,\n      \
solpar,reffor,conpar,state,options,soldim]];\n   If [neval==3,\n     \
{nextsol,status}=IncStep3[ptags,geopar,matpar,fabpar,\n      \
solpar,reffor,conpar,state,options,soldim]];\n   If [neval==4,\n     \
{nextsol,status}=IncStep4[ptags,geopar,matpar,fabpar,\n      \
solpar,reffor,conpar,state,options,soldim]];\n   Return[{nextsol,status}]]; "
}], "Input",
  CellFrame->True,
  CellMargins->{{14, 37}, {Inherited, Inherited}},
  InitializationCell->True,
  CellLabelMargins->{{8, Inherited}, {Inherited, Inherited}},
  ImageRegion->{{-0, 1}, {0, 1}},
  Background->RGBColor[0, 1, 0]],

Cell[TextData[{
  StyleBox["SingularityFix",
    FontWeight->"Bold"],
  ": Tries to recover from a singular stiffness matrix by perturbing the \
state \nwith small random numbers. Gives up after 3 unsuccesful tries"
}], "Text",
  CellFrame->True,
  CellMargins->{{14, 37}, {Inherited, Inherited}},
  CellLabelMargins->{{8, Inherited}, {Inherited, Inherited}},
  ImageRegion->{{-0, 1}, {0, 1}},
  Background->RGBColor[1, 1, 0]],

Cell[TextData[{
  "SingularityFix",
  StyleBox["[",
    FontColor->GrayLevel[0]],
  "ptags_",
  StyleBox[",geopar_,matpar_,fabpar_,solpar_,",
    FontColor->GrayLevel[0]],
  "reffor_,\n    conpar_,state_,options_,soldim_",
  StyleBox["]:=Module[",
    FontColor->GrayLevel[0]],
  "{nextsol,q,nudge,\n    n,m,un,OK,ueps=10.^(-10),numer=True},\n    nudge=0; \
SeedRandom[7654321]; ",
  StyleBox["{n,m}=soldim;",
    FontColor->GrayLevel[0]],
  "\n    nextsol=Table[\"Null\",{m}]; OK=False;\n    For [nudge=1, nudge<=3, \
nudge++, \n        K=TanStiffOfMisesTruss[sprop,un]; K=N[K];",
  StyleBox["\n",
    FontColor->RGBColor[1, 0, 0]],
  "       {nextsol,status}=IncStep[ptags,geopar,matpar,fabpar,\n        \
solpar,reffor,conpar,state,options",
  ",soldim]",
  "; \n        If [StringTake[status,8]!=\"Singular\",\n        OK=True; \
Break[]];\n        un+=Table[Random[]*ueps,{2}] \n        ];\n    If [OK, \
Return [{N[v],N[q],\" \"}]];\n    Return [{Null,\"Cant escape \
singularity\"}]];"
}], "Input",
  CellFrame->True,
  CellMargins->{{14, 37}, {Inherited, Inherited}},
  InitializationCell->True,
  CellLabelMargins->{{8, Inherited}, {Inherited, Inherited}},
  ImageRegion->{{-0, 1}, {0, 1}},
  Background->RGBColor[0, 1, 0]],

Cell[TextData[{
  StyleBox["NonLinSolver: ",
    FontWeight->"Bold",
    FontColor->GrayLevel[0]],
  StyleBox["Solution advancing driver: checks data; if OK sets ref state & \
runs analysis",
    FontColor->GrayLevel[0]]
}], "Text",
  CellFrame->True,
  CellMargins->{{14, 37}, {Inherited, Inherited}},
  CellLabelMargins->{{8, Inherited}, {Inherited, Inherited}},
  ImageRegion->{{-0, 1}, {0, 1}},
  Background->RGBColor[1, 1, 0]],

Cell[TextData[{
  "NonLinSolver[ptags_,geopar_,matpar_,fabpar_,solpar_,reffor_,\n    \
options_]:=Module[{conpar,state,problem,int,ics, \n     \
n,m,nmax,\[Lambda]max,umax,ell,ellmin,ellmax,acctol,cscale,A,c,w,\n     \
soltab,solnull,numdof,icsOK,\n     \
\[Lambda]0,u0,v0,q0,r0,Kdet0,sol0,status=\" \"},\n   \
{nmax,\[Lambda]max,umax,ell,ellmin,ellmax,acctol,cscale}=solpar;\n    \
solnull={0,0,0,0,0,0,0,0,0,0,\" \"};\n   {problem,int,ics}=ptags; \
numdof=ProblemInfo[problem]; \n    If [numdof<=0, status=\"Unimplemented \
problem\";\n                   Return[{solnull,status}]];\n   \
{A,c,w}=ButcherTable[int]; \n    If [c==Null,   status=\"Unimplemented int\"; \
\n                   Return[{solnull,status}]];\n    \
icsOK=(ics==\"LC\")||(ics==\"SC\")||(ics==\"AC\")||(ics==\"HC\");  \n    If \
[!icsOK,    status=\"Unimplemented ics\"; \n                   \
Return[{solnull,status}]];\n    n=0; \[Lambda]0=0; conpar={\[Lambda]0}; \
u0=Table[0,{numdof}]; state=u0; \n    {{v0,q0,Kdet0},status}= \
IncVel[ptags,geopar,\n         \
matpar,fabpar,solpar,reffor,conpar,state,options];\n    If [status!=\" \", \
Return[{solnull,status}]];\n    r0=TotalResidual[ptags,geopar,matpar,fabpar,\n\
         solpar,reffor,conpar,state,options];",
  StyleBox[" ",
    FontColor->RGBColor[1, 0, 0]],
  "\n    soltab=Table[solnull,{nmax+2}]; \n    \
sol0={n,\[Lambda]0,u0,v0,r0,q0,Kdet0,ell,1,1,\"Ref config\"};\n    \
soltab[[1]]=sol0; m=Length[sol0]; \n    While [status==\" \" && n<=nmax, \n   \
   {nextsol,status}=IncStep[ptags,geopar,matpar,\n       \
fabpar,solpar,reffor,conpar,state,options,{n,m}]; \n       n++; \
conpar={nextsol[[2]]}; state=nextsol[[3]]; \n       soltab[[n+1]]=nextsol;",
  StyleBox[" ",
    FontColor->RGBColor[1, 0, 0]],
  "\n          ];\n    Return[{Take[soltab,n],status}]]; "
}], "Input",
  CellFrame->True,
  CellMargins->{{14, 37}, {Inherited, Inherited}},
  InitializationCell->True,
  CellLabelMargins->{{8, Inherited}, {Inherited, Inherited}},
  ImageRegion->{{-0, 1}, {0, 1}},
  Background->RGBColor[0, 1, 0]],

Cell[TextData[StyleBox["Benchmark problem scripts",
  FontWeight->"Bold",
  FontColor->GrayLevel[0]]], "Text",
  CellFrame->True,
  CellMargins->{{14, 37}, {Inherited, Inherited}},
  CellLabelMargins->{{8, Inherited}, {Inherited, Inherited}},
  ImageRegion->{{-0, 1}, {0, 1}},
  FontColor->RGBColor[0, 1, 1],
  Background->RGBColor[0, 1, 1]],

Cell["\<\
Driver script to define and run the Mises Truss (MT) with:   S=2, \
H=1, Em=1, A0=1, base load
q={0,-Em*A0},  fixed steplength ell usually .1 to .001.  Continuation scheme: \
defined in ptags.\
\>", "Text",
  CellFrame->True,
  CellMargins->{{14, 37}, {Inherited, Inherited}},
  CellLabelMargins->{{8, Inherited}, {Inherited, Inherited}},
  ImageRegion->{{-0, 1}, {0, 1}},
  Background->RGBColor[1, 1, 0]],

Cell[CellGroupData[{

Cell["\<\
ClearAll[S,H,Em,A0,nmax,ell,cscale];
 
(* Set tags: {benchmark-problem,advancing-integration-scheme,
   increment-control-strategy}. Internally: problem,int,ics.
   Note that ics links w/ cscale, declared as {} for default *)

ptags={\"MT\",\"FE\",\"LC\"}; cscale={};

(* Define structural properties: {span height H,
    modulus Em, cross section A0} *)

S=2; H=1; Em=1; A0=1; 
geopar={S,H}; matpar={Em}; fabpar={A0,A0,0};

(* Ext force is \[Lambda] \[Times] reference force: \
\[Lambda]*reffor=\[Lambda]*{0,-Em*A0} *)

reffor={0,-Em*A0};

(* Define solution parameters - see above re cscale  *)  

nmax=40; \[Lambda]max=2; umax={2*S,2*H}; ell=0.04; acctol=0; cscale={}; 
solpar={nmax,\[Lambda]max,umax,ell,ell,ell,acctol,cscale};

(*  Specify numerical (floating-point) computation *)

numer=True; options={numer};

(*  Carry out nonlinear solution, which returns in soltab *)

{soltab,status}=NonLinSolver[ptags,geopar,
    matpar,fabpar,solpar,reffor,options];

(*  Print solution table *)

doflab={\"uX\",\"uY\"}; digits={{5,3}};
PrintPISolTable[ptags,soltab,doflab,digits];
If [status!=\" \",Print[status]];

(*  Response plots *)

imgsiz=350;  (* plot frame width is 350 pt *)
PlotPISolTable[ptags,soltab,{\"uY\",\"\[Lambda]\"},doflab,imgsiz];
PlotPISolTable[ptags,soltab,{\"||r||\",\"\[Lambda]\"},doflab,imgsiz]; 
PlotPISolTable[ptags,soltab,{\"step\",\"||r||\"},doflab,imgsiz];  
\
\>", "Input",
  CellFrame->True,
  CellMargins->{{14, 37}, {Inherited, Inherited}},
  CellLabelMargins->{{8, Inherited}, {Inherited, Inherited}},
  ImageRegion->{{-0, 1}, {0, 1}},
  Background->RGBColor[0, 1, 0]],

Cell["\<\
Changed int from FE to EMR, and ics from LC to AC - things should \
not look so bad\
\>", "Text",
  CellFrame->True,
  CellMargins->{{14, 37}, {Inherited, Inherited}},
  CellLabelMargins->{{8, Inherited}, {Inherited, Inherited}},
  ImageRegion->{{-0, 1}, {0, 1}},
  Background->RGBColor[1, 1, 0]],

Cell["\<\
ClearAll[S,H,Em,A0,nmax,ell,cscale];
 
(* Set tags: {benchmark-problem,advancing-integration-scheme,
   increment-control-strategy}. Internally: problem,int,ics.
   Note that ics links w/ cscale, declared as {} for default *)

ptags={\"MT\",\"EMR\",\"AC\"}; cscale={};

(* Define structural properties: {span height H,
    modulus Em, cross section A0} *)

S=2; H=1; Em=1; A0=1; 
geopar={S,H}; matpar={Em}; fabpar={A0,A0,0};

(* Ext force is \[Lambda] \[Times] reference force: \
\[Lambda]*reffor=\[Lambda]*{0,-Em*A0} *)

reffor={0,-Em*A0};

(* Define solution parameters - see above re cscale  *)  

nmax=60; \[Lambda]max=2; umax={2*S,2*H}; ell=0.04; acctol=0; cscale={}; 
solpar={nmax,\[Lambda]max,umax,ell,ell,ell,acctol,cscale};

(*  Specify numerical (floating-point) computation *)

numer=True; options={numer};

(*  Carry out nonlinear solution, which returns in soltab *)

{soltab,status}=NonLinSolver[ptags,geopar,
    matpar,fabpar,solpar,reffor,options];

(*  Print solution table *)

doflab={\"uX\",\"uY\"}; digits={{5,3}};
PrintPISolTable[ptags,soltab,doflab,digits];
If [status!=\" \",Print[status]];

(*  Response plots *)

imgsiz=350;  (* plot frame width is 350 pt *)
PlotPISolTable[ptags,soltab,{\"uY\",\"\[Lambda]\"},doflab,imgsiz];
PlotPISolTable[ptags,soltab,{\"||r||\",\"\[Lambda]\"},doflab,imgsiz]; 
PlotPISolTable[ptags,soltab,{\"step\",\"||r||\"},doflab,imgsiz];  

\
\>", "Input",
  CellFrame->True,
  CellMargins->{{14, 37}, {Inherited, Inherited}},
  CellLabelMargins->{{8, Inherited}, {Inherited, Inherited}},
  ImageRegion->{{-0, 1}, {0, 1}},
  Background->RGBColor[0, 1, 0]]
}, Open  ]]
},
FrontEndVersion->"4.2 for Macintosh",
ScreenRectangle->{{0, 1920}, {0, 1180}},
AutoGeneratedPackage->None,
WindowToolbars->{},
CellGrouping->Manual,
WindowSize->{1466, 1083},
WindowMargins->{{97, Automatic}, {Automatic, 7}},
PrivateNotebookOptions->{"ColorPalette"->{RGBColor, -1}},
ShowCellLabel->True,
ShowCellTags->False,
RenderingOptions->{"ObjectDithering"->True,
"RasterDithering"->False},
Magnification->1.5,
MacintoshSystemPageSetup->"\<\
00<0001804P000000]P2:?oQon82n@960dL5:0?l0080001804P000000]P2:001
0000I00000400`<300000BL?00400@00000000000000060801T1T00000000000
00000000000000000000000000000000\>"
]

(*******************************************************************
Cached data follows.  If you edit this Notebook file directly, not
using Mathematica, you must remove the line containing CacheID at
the top of  the file.  The cache data will then be recreated when
you save this file from within Mathematica.
*******************************************************************)

(*CellTagsOutline
CellTagsIndex->{}
*)

(*CellTagsIndex
CellTagsIndex->{}
*)

(*NotebookFileOutline
Notebook[{
Cell[1754, 51, 352, 7, 138, "Section"],
Cell[2109, 60, 3391, 76, 913, "Text"],
Cell[5503, 138, 307, 7, 72, "Text"],
Cell[5813, 147, 349, 9, 71, "Text"],
Cell[6165, 158, 414, 11, 122, "Input",
  InitializationCell->True],
Cell[6582, 171, 617, 16, 96, "Text"],
Cell[7202, 189, 1589, 28, 682, "Input",
  InitializationCell->True],
Cell[8794, 219, 382, 11, 72, "Text"],
Cell[9179, 232, 469, 12, 122, "Input",
  InitializationCell->True],
Cell[9651, 246, 606, 16, 123, "Text"],
Cell[10260, 264, 2072, 52, 662, "Input",
  InitializationCell->True],
Cell[12335, 318, 333, 9, 72, "Text"],
Cell[12671, 329, 633, 15, 146, "Text"],
Cell[13307, 346, 3641, 59, 1202, "Input",
  InitializationCell->True],
Cell[16951, 407, 341, 9, 72, "Text"],
Cell[17295, 418, 3198, 58, 1042, "Input",
  InitializationCell->True],
Cell[20496, 478, 398, 10, 72, "Text"],
Cell[20897, 490, 625, 13, 171, "Text"],
Cell[21525, 505, 2412, 42, 1062, "Input",
  InitializationCell->True],
Cell[23940, 549, 650, 13, 171, "Text"],
Cell[24593, 564, 1722, 31, 602, "Input",
  InitializationCell->True],
Cell[26318, 597, 342, 9, 72, "Text"],
Cell[26663, 608, 383, 11, 122, "Input",
  InitializationCell->True],
Cell[27049, 621, 398, 10, 72, "Text"],
Cell[27450, 633, 342, 9, 72, "Text"],
Cell[27795, 644, 2172, 42, 622, "Input",
  InitializationCell->True],
Cell[29970, 688, 330, 9, 72, "Text"],
Cell[30303, 699, 2223, 39, 702, "Input",
  InitializationCell->True],
Cell[32529, 740, 335, 9, 72, "Text"],
Cell[32867, 751, 2155, 42, 622, "Input",
  InitializationCell->True],
Cell[35025, 795, 395, 10, 72, "Text"],
Cell[35423, 807, 491, 13, 96, "Text"],
Cell[35917, 822, 3988, 79, 1102, "Input",
  InitializationCell->True],
Cell[39908, 903, 628, 17, 99, "Text"],
Cell[40539, 922, 4499, 112, 942, "Input",
  InitializationCell->True],
Cell[45041, 1036, 762, 19, 123, "Text"],
Cell[45806, 1057, 3046, 87, 582, "Input",
  InitializationCell->True],
Cell[48855, 1146, 751, 19, 123, "Text"],
Cell[49609, 1167, 3721, 102, 702, "Input",
  InitializationCell->True],
Cell[53333, 1271, 753, 19, 123, "Text"],
Cell[54089, 1292, 4407, 119, 842, "Input",
  InitializationCell->True],
Cell[58499, 1413, 590, 19, 99, "Text"],
Cell[59092, 1434, 1443, 35, 422, "Input",
  InitializationCell->True],
Cell[60538, 1471, 426, 10, 96, "Text"],
Cell[60967, 1483, 1232, 31, 342, "Input",
  InitializationCell->True],
Cell[62202, 1516, 431, 12, 72, "Text"],
Cell[62636, 1530, 2044, 37, 682, "Input",
  InitializationCell->True],
Cell[64683, 1569, 341, 8, 72, "Text"],
Cell[65027, 1579, 415, 10, 93, "Text"],

Cell[CellGroupData[{
Cell[65467, 1593, 1613, 52, 942, "Input"],
Cell[67083, 1647, 306, 8, 69, "Text"],
Cell[67392, 1657, 1615, 53, 962, "Input"]
}, Open  ]]
}
]
*)



(*******************************************************************
End of Mathematica Notebook file.
*******************************************************************)