NullNullNullNull

Off[General::"spell"]; Off[General::"spell1"]; $Line = 0;

BeginPackage["MathEcon`"]; Needs["Graphics`Arrow`"];

Unprotect[fibprime, directed, vector2d, directed3d, vector3d, 
      addition2d, subtraction2d, projection2d, projection3d, rotation2d, 
      rank, rowop, colop, submatrix, minor, cofactor, aug, colrep, singQ, 
      signQL, borderB, gradf, hessian, conhess, arc2d, tangentline, fjohn, 
      separableode, linearode, tangentfield, wronskian, wis, trajectory]; 
    Clear[fibprime, directed, vector2d, directed3d, vector3d, addition2d, 
      subtraction2d, projection2d, projection3d, rotation2d, rank, rowop, 
      colop, submatrix, minor, cofactor, aug, colrep, singQ, signQL, borderB, 
      gradf, hessian, conhess, arc2d, tangentline, fjohn, separableode, 
      linearode, tangentfield, wronskian, wis, trajectory];

MathEcon::"usage" = 
      "MathEcon is a package of functions that is used with the book:\r\r \
Mathematics and Mathematica for Economists \r         by Cliff Huang and \
Philip Crooke.\r\r It contains the following functions: \r\r fibprime  \
directed  directed \r vector2d vector3d  addition2d \r subtraction2d  \
projection2d \r rotation2d  rank  rowop \r  colop  submatrix  minor \r \
cofactor  aug  colrep \r singQ  signQL  borderB \r gradf  hessian  conhess \r \
arc2d  tangentline  fjohn \r separableode  linearode tangentfield \r \
wronskian  wis  trajectory \r";

fibprime::"usage" = 
      "fibprime[n] calculates the prime Fibonnaci numbers that are less \
than or equal to an integer n.  For example, \r\rfibprime[37].";

directed::"usage" = 
      "directed[pt1,pt2,\"label\",opt] generates a 2D Graphics object for a \
directed line segment that begins at the point pt1 and ends at the point pt2. \
It allows the line segment to have a \"label\" (a string), and a dashed line \
style is used when the fourth argument is present. For example, \r\r\
directed[{-1,2},{3,2},\"a\",dash].";

vector2d::"usage" = 
      "vector2d[pt,\"label\",opt] generates a 2D Grapbics object for a \
vector from the origin to the point specified in the first argument. It \
allows the vector to have a label (a string) and a dashed line style is used \
when the third argument is present.  For example, \r\r\
vector[{3,4},\"a\",dash].";

directed3d::"usage" = 
      "directed3d[pt1,pt2,\"label\",opt] generates a 3D Graphics object for \
a directed line segment that begins at the point pt1 and ends at the point \
pt2. It allows the line segment to have a label (a string) and a dashed line \
style can be used when the fourth argument is present.  For example, \r\r\
directed3d[{0,1,-1},{3,4,-2},\"a\",dash].";

vector3d::"usage" = 
      "vector3d[pt,\"label\",opt] generates a 3D Grapbics object for a \
vector from the origin to the point specified in the first argument. It \
allows the vector to have a label (string) and a dashed line style is used \
when the third argument is present.  For example, \r\r\
vector3d[{1,-9,0},\"a\",dash].";

addition2d::"usage" = 
      "addition2d[pt1,\"label1\",pt2,\"label2\"] generates a 2D Graphics \
object that illustrates the addition of two vectors from the origin that \
terminate at the points pt1 and pt2, respectively.  Each vector has a label \
(string) and the vector that is the sum is labeled.  For example, \r\r\
addition2d[{-1,2},\"a\",{0,-1},\"b\"].";

subtraction2d::"usage" = 
      "subrtraction2d[pt1,\"label1\",pt2,\"label2\"] generates a 2D \
Graphics object that illustrates the difference of two vectors from the \
origin that terminate at the points pt1 and pt2, respectively.  Each vector \
has a label (a string) and the vector that is the difference is labeled.  For \
example, \r\rsubtraction2d[{-1,2},\"a\",{0,-1},\"b\"].";

projection2d::"usage" = 
      "projection2d[pt1,\"label1\",pt2] generates a 2D Graphics object that \
illustrates the projection of a vector that starts at the origin and ends at \
pt1 onto another vector that starts at the origin and ends at pt2.  The first \
vector can have a label (a string).  For example, \r\rprojection2d[{-2,1},\"a\
\",{1,4}].";

projection3d::"usage" = 
      "projection3d[pt1,\"label1\",pt2] generates a 3D Graphics object that \
illustrates the projection of a vector that starts at the origin and ends at \
pt1 onto another vector that starts at the origin and ends at pt2. The first \
vector can have a label (a string).  For example, \r\r\
projection3d[{1,0,2},\"a\",{3,1,2}].";

rotation2d::"usage" = 
      "rotation2d[pt,\"label\",ang] generates a 2D Graphics object that \
illustrates the rotation of the vector that starts at the origin and ends at \
the point pt through the angle ang.  The vector can have a label (a string).  \
For example, \r\rrotation2d[{-1,2},\"a\",Pi/4]."; rank::"usage" = 
      "rank[lst] calculates the dimension of the space spanned by the \
vectors contained in its argument (list).  For example, \r\r\
rank[{{1,-1,9},{1,-2,-3}}].";

rowop::"usage" = 
      "rowop[rowlist] generates the elementary row operation matrix that \
are dictated by its argument which is a list of elemntary row operations that \
are to be performed.  For example, \r\rrowop[{row[3],row[2],row[1]}].";

colop::"usage" = 
      "colop[collist] generates the elementary column operation matrix that \
are dictated by its argument which is a list of elemntary column operations \
that are to be performed.  For example, \r\r \
colop[{col[2],col[3],col[1]}]."; row::"usage" = 
      "row[i] is a function used in the rowop function to represent the ith \
row of the identity matrix.  For example, \r\rrow[3]."; col::"usage" = 
      "col[i] is a function used in the colop function to represent the ith \
column of the identity matrix.  For example,  \r\rcol[2].";

submatrix::"usage" = 
      "submatrix[mat,irow,jcol] calculates the submatrix of a given matrix \
mat by removing row irow and column jcol.  For example, \r\r\
submatrix[{{1,9},{6,-5}},1,2]."; minor::"usage" = 
      "minor[mat,irow,jcol] calculates the (irow,jcol) minor of a given \
matrix mat.  For example, \r\rminor[{{1,9},{6,-5}},1,2].";

cofactor::"usage" = 
      "cofactor[mat,irow,jcol] calculates the (irow,jcol) cofactor of a \
given matrix mat.  For example, \r\rcofactor[{{1,9},{6,-5}},1,2].";

aug::"usage" = 
      "aug[mat,vect] forms the augmented matrix of the matrix mat and the \
vector vect.  For example, \r\raug[{{1,2},{0,1}},{-3,2}]."; colrep::"usage" = 
      "colrep[mat,vect,icol] replaces the icol column of a given matrix mat \
by the vector vect.  For example, \r\rcolrep[{{1,2},{0,1}},{-3,2},2].";

changebasis::"usage" = 
      "changebasis[vect1,bas,\"label1\",\"basislabel\"] creates a graphic \
that illustrates how the vector vect1 given in the usual basis {{1,0},{0,1}} \
appears in the basis bas.  The given vector and the basis vectors can have \
labels (strings).  For example, \r\r\
changebasis[{-1,3},{{-1,1},{0,1}},\"a\",\"b\"].";

signQ::"usage" = 
      "signQ[mat] creates a list of the principal minors of the matrix mat. \
 For example, \r\rsignQ[{{2,3,1},{4,3,1},{0,9,-3}}].";

signQL::"usage" = 
      "signQL[mat] creates a list of the leading principal minors of the \
matrix mat.  For example, \r\rsignQL[{{2,3,1},{4,3,1},{0,9,-3}}].";

borderB::"usage" = 
      "borderB[mat,vect] genrates the border matrix of the matrix mat and \
the vector vect.  For example, \r\r \
borderB[{{2,3,1},{4,3,1},{0,9,-3}},{-1,0,4}].";

gradf::"usage" = 
      "gradf[fct,lst] calculates the gradient of the function fct with \
respect to the variables specified by lst (a list).  For example, \r\r\
gradf[x*y-y^2*z+w,{x,y,z,w}].";

hessian::"usage" = 
      "hessian[fct,lst] generates the hessian matrix of the function fct \
with respect to the variables specified by lst (a list).  For example, \r\r\
hessian[x*y*w-y^2*z+w^3,{x,y,z,w}].";

conhess::"usage" = 
      "conhess[fct,lst1,fctlst,lst2] creates the second-order differential \
(a matrix) of the function fct with respect to the variables specified by \
lst1 (a list) where the variables in the list lst2 are constrained by \
functions specified in the list of functions fctlst.  For example, \r\r\
conhess[x*y*w-y^2*z+w^3,{x,y,z,w},{x+y*z},{w}].";

arc2d::"usage" = 
      "arc2d[fct,dom,seg,opt] generates a 2D Graphics object that \
illustrates the graph of the function fct over the interval dom (a list) with \
the graph highlighted on the subinterval seg (a list).  The fourth argument \
is optional.  If the fourth argument is chord, then a chord is drawned \
between points specified by seg. For example, \r\r\
arc2d[x^2-1,{x,-3,4},{x,0,1},chord].";

tangentline::"usage" = 
      "tangentline[fct,int,pt] generates a 2D Graphics object that contains \
the graph of the function fct over the interval int (a list) with the tangent \
line to the graph of fct at the point pt. The graphics can then be display \
with the Show function. For example, \r\r \
Show[tangentline[Cos[2*Pi*x],{x,-1,1},Pi/4],Axes->True];"; fjohn::"usage" = 
      "fjohn[pt,const,objfct,var,options] generates a 2D Graphics to \
illustrate the gradients of an objective function objfct and binding \
constraint const (a list) at the point pt in a nonlinear programming problem. \
 The fourth argument of fjohn contains a list of the variables for the \
objective and constraint functions.  The last argument is optional.  If it is \
the string \"label\", then the gradient vectors are labeled; if it is the \
string \"cone\", then the gradients of the binding contraints and the \
negative of the gradinet of the objective function is drawn and labeled. For \
example, \r\rfjohn[{1,1},{x-y,x^2-y},x^2+x^3,{x,y},\"label\"]."; separableode::"usage" = 
      "separableode[var1,var2,fct1,fct2] attempts to find the general \
implicit and explicit solutions of the separable, first order differential \
equation: fct1 d(var1) + fct2 d(var2) = 0. The solutions are presented as a \
list of rules.  The arbitrary constant is denoted by C[1].  For example, \r\r\
separableode[x,y,Cos[x],-y*Sin[y^2]]."; linearode::"usage" = 
      "linearode[var1,fct1,fct2] attempts to find the general, explicit \
solution of the first order, linear ODE: d(var2)/d(var1) + fct1(var1) var2 = \
fct2(var1).  The output is an expression involving the symbol var1 and the \
arbitrary constant C[1].  For example, \r\rlinearode[x,1/x,x^2]."; tangentfield::"usage" = 
      "tangentfield[fct,var1range,var2range,int,opt] generates a 2D \
Graphics object that illustrate the tangent field of the differential \
equation: d(var2)/d(var1) = fct where fct is a function of var1 and var2, \
var1range is a list that is the range of the variable var1, and var2range is \
a list for the range of the variable var2 over which the tangent field will \
be displayed.  The fourth argument int is a positive integer that is related \
to the number of tangent lines to be drawn. The last argument is optional and \
can include 2D Graphics options.  For example, \r\r\
tangentfield[x*y^3,{x,-1,2},{y,0,2},10,AxesLabel->{\"x\",\"y\"}]."; wronskian::"usage" = 
      "wronskian[fctlst,var] computes the Wronkian determinant of a list of \
functions in the first argument whose independent variable is the second \
argument var.  For example, \r\rwronskian[{1,x,x^2},x]."; wis::"usage" = 
      "wis[fct,fctlst,var] computes the determinants that are used in the \
Variation of Parameters Method. For example, \r\rwis[Exp[x],{1,x,x^2},x]."; 
    trajectory::"usage" = 
      "trajectory[fcts,dom,opt] displays the trajectory of the list of \
functions fcts given on the domain dom (a list).  The last argument is \
optional and can contain any option for ParametricPlot. Arrows are placed on \
the graph of the parametric functions to indicate the direction of increasing \
values of the independent variable.  For example, \r\r\
trajectory[{x^2+1,Cos[x]},{x,-1,3}, AxesLabel->{\"y1\",\"y2\"}].";

Begin["`Private`"]; fibprime[n_Integer] := 
      Module[{k, t, lessprime, x}, Clear[f, m]; f[0] = 1; f[1] = 1; 
        f[m_ /; m ≥ 2] := f[m] = f[m - 1] + f[m - 2]; 
        t = Table[f[k], {k, 0, n}]; 
        lessprime[x_] := PrimeQ[x] && x ≤ n; 
        Select[t, lessprime]]; directed[start_List, end_List, label_String, dash___] := 
      Module[{p, f, aa, bb, cc, dd, len, setdash, A, B, temp, end1, end2}, {
          p = (π/10); f = 0.12; setdash = {N[(1/72)]}; 
          If[dash == Null, setdash = {}]; 
          aa = start[[1]]; 
          bb = start[[2]]; 
          cc = end[[1]]; 
          dd = end[[2]]; 
          len = (Sqrt[(((aa - cc))^2) + (((bb - dd))^2)]); 
          A = {{Cos[p], Sin[p]}, {-Sin[p], Cos[p]}}; 
          B = {{Cos[-p], Sin[-p]}, {-Sin[-p], Cos[-p]}}; 
          temp = end - start; end1 = end - ((f A . temp)/len); 
          end2 = end - ((f B . temp)/len); 
          Graphics[{AbsoluteThickness[N[(1/72)]], AbsoluteDashing[{}], 
              Line[{end, end1}], Line[{end, end2}], Dashing[setdash], 
              Line[{start, end}], 
              Text[FontForm[label, {"Courier-Bold", 12}], {
                  ((aa + cc)/2), ((bb + dd)/2)}, {0, 2}]}]}]; vector2d[v_List, lab_String, styl___] := 
      directed[{0, 0}, v, lab, styl]; directed3d[start_List, end_List, label_String, dash___] := 
      Module[{p, hd, len, AA, BB, aa, bb, cc, dd, ee, ff, setdash, end1, 
          end2}, {p = (π/10); hd = 0.25; setdash = {N[(1/72)]}; 
          If[dash == Null, setdash = {}]; 
          aa = start[[1]]; 
          bb = start[[2]]; 
          cc = start[[3]]; 
          dd = end[[1]]; 
          ee = end[[2]]; 
          ff = end[[3]]; 
          len = (Sqrt[(((aa - dd))^2) + (((bb - ee))^2) + (((cc - ff))^2)]); 
          AA = {{Cos[p], Sin[p], 0}, {-Sin[p], Cos[p], 0}, {0, 0, 1}}; 
          BB = {{Cos[-p], Sin[-p], 0}, {-Sin[-p], 
                Cos[-p], 0}, {0, 0, 1}}; 
          end1 = end - ((hd AA . (end - start))/len); 
          end2 = end - ((hd BB . (end - start))/len); 
          Graphics3D[{AbsoluteThickness[N[(9/72)]], Dashing[{}], 
              Line[{end, end1}], Line[{end, end2}], Dashing[setdash], 
              Line[{start, end}], 
              Text[FontForm[label, {"Geneva-Bold", 12}], {((aa + dd)/2), 
                  ((bb + ee)/2), ((cc + ff)/2)}, {0, 2}]}]}]vector3d[v_List, label_String, dash___] := 
      directed3d[{0, 0, 0}, v, label, dash]; addition2d[a_List, labela_String, b_List, labelb_String] := 
      Module[{sumab}, {directed[{0, 0}, a, labela], 
          directed[{0, 0}, b, labelb], directed[a, b + a, " ", dash], 
          directed[{0, 0}, a + b, StringJoin[{labela, "+", labelb}]]}]; 
    subtraction2d[a_List, labela_String, b_List, labelb_String] := 
      Module[{subab}, {directed[{0, 0}, a, labela], 
          directed[{0, 0}, b, labelb], directed[a, a - b, " ", dash], 
          directed[{0, 0}, a - b, StringJoin[{labela, "-", labelb}]]}]; 
    projection2d[a_List, labela_String, b_List] := 
      Module[{temp}, {directed[{0, 0}, a, labela], 
          directed[{0, 0}, b, ""], 
          directed[{0, 0}, ((a . b b)/b) . b, "proj"], 
          Graphics[{Dashing[{0.01}], Line[{a, ((a . b b)/b) . b}]}]}]; 
    projection3d[a_List, labela_String, b_List] := 
      Module[{temp}, {directed3d[{0, 0, 0}, a, labela], 
          directed3d[{0, 0, 0}, b, ""], 
          directed3d[{0, 0, 0}, ((a . b b)/b) . b, "proj"], 
          Graphics3D[{Dashing[{0.02}], Line[{a, ((a . b b)/b) . b}]}]}]; 
    rotation2d[a_List, labela_String, angle_] := 
      Module[{A, temp, p}, {
          A = {{Cos[angle], -Sin[angle]}, {Sin[angle], Cos[angle]}}; 
          temp = A . a; 
          p = ArcTan[a[[1]], 
              a[[2]]]; {
            directed[{0, 0}, a, labela], directed[{0, 0}, temp, "rot"], 
            Graphics[{Dashing[{0.02}], 
                Circle[{0, 0}, (Sqrt[a . a]), {p, p + angle}]}]}}]; rank[v_List] := Length[Transpose[v]] - Length[NullSpace[v]]; Clear[rowop]; rowop[x_List] := 
      Module[{n}, n = Length[x]; 
        row[y_] := 
          IdentityMatrix[n][[y]]; 
        x]; colop[x_List] := 
      Module[{n}, n = Length[x]; 
        col[y_] := 
          IdentityMatrix[n][[y]]; 
        Transpose[x]]; submatrix[A_List, irow_Integer, jcol_Integer] := 
      Module[{index}, index = Table[k, {k, Length[A]}]; 
        A[[Drop[index, {irow, irow}], 
          Drop[index, {jcol, jcol}]]]]; minor[A_List, irow_Integer, jcol_Integer] := 
      Module[{index}, index = Table[k, {k, Length[A]}]; 
        Det[A[[Drop[index, {irow, irow}], 
            Drop[index, {jcol, jcol}]]]]]; cofactor[A_List, irow_Integer, jcol_Integer] := 
      Module[{index}, index = Table[k, {k, Length[A]}]; 
        (((-1))^(irow + jcol)) 
          Det[A[[Drop[index, {irow, irow}], 
              Drop[index, {jcol, jcol}]]]]]; aug[a_List, v_List] := Transpose[Append[Transpose[a], v]]; colrep[a_List, v_List, i_Integer] := 
      Transpose[ReplacePart[Transpose[a], v, i]]; Clear[changebasis]; changebasis[x_List, p_List, xlabel_String: "x", 
        plabel_String: "p"] := 
      Module[{origin, y, points, lines, vectx, vectp1, vectp2, vectI1, 
          vecI2}, origin = {0, 0}; y = Inverse[p] . x; 
        points = 
          Graphics[{PointSize[0.02], Point[x], Point[{1, 0}], Point[{0, 1}], 
              Point[{x[[1]], 0}], 
              Text[xlabel<>"1", {
                  x[[1]], 0}, {1, 
                  -1}], 
              Point[{0, x[[2]]}], 
              Text[xlabel<>"2", {0, 
                  x[[2]]}, {-1, 1}], 
              Point[p[[1]]], 
              Point[p[[2]]], 
              Point[y[[1]] 
                  p[[1]]], 
              Text[" y1", 
                y[[1]] 
                  p[[1]], {-1, 0}], 
              Point[y[[2]] 
                  p[[2]]], 
              Text["y2 ", 
                y[[2]] 
                  p[[2]], {1, 0}]}]; 
        lines = Graphics[{
              Line[{origin, 
                  y[[1]] 
                    p[[1]]}], 
              Line[{origin, 
                  y[[2]] 
                    p[[2]]}], 
              Dashing[{0.01}], 
              Line[{x, {x[[1]], 0}}], 
              Line[{x, {0, x[[2]]}}], 
              Line[{x, 
                  y[[1]] 
                    p[[1]]}], 
              Line[{x, 
                  y[[2]] 
                    p[[2]]}]}]; 
        vectx = vector2d[x, xlabel]; 
        vectp1 = 
          vector2d[p[[1]], 
            plabel<>"1"]; 
        vectp2 = 
          vector2d[p[[2]], 
            plabel<>"2"]; vectI1 = vector2d[{1, 0}, ""]; 
        vectI2 = vector2d[{0, 1}, ""]; {points, lines, vectx, vectp1, 
          vectp2, vectI1, vectI2}]; signQ[A_List] := 
      Module[{k}, k = Length[A]; 
        Table[Minors[A, i][[j, j]], {
            i, k}, {j, Binomial[k, i]}]]; signQL[A_List] := 
      Module[{k}, k = Length[A]; 
        Table[Minors[A, i][[1, 1]], {
            i, k}]]; borderB[A_List, b_List] := 
      Module[{cons, consrow, conscol}, cons = If[VectorQ[b], {b}, b]; 
        consrow = 0 IdentityMatrix[Length[cons]]; 
        Do[consrow[[i]] = 
            Join[consrow[[i]], 
              cons[[i]]], {i, 
            Length[cons]}]; conscol = Transpose[cons]; 
        Do[conscol[[i]] = 
            Join[conscol[[i]], 
              A[[i]]], {i, Length[A]}]; 
        bord = Join[consrow, conscol]; bord]; Clear[gradf]; gradf[f_, v_List] := 
      Module[{n, i}, n = Length[v]; 
        Table[(∂_(v[[i]]))f, {
            i, 1, n}]]; hessian[f_, v_List] := 
      Module[{n, i, j}, n = Length[v]; 
        Table[(∂_(v[[i]], 
              v[[j]]))f, {i, 1, n}, {j, 
            1, n}]]; conhess[f_, v_List, g_, u_List] := 
      Module[{m, ch, k}, m = Length[u]; ch = hessian[f, v]; 
        Do[ch = ch + 
              (∂_(u[[k]]))f 
                hessian[g[[k]], v], {k, 
            1, m}]; ch]; signQL[A_List] := 
      Module[{k}, k = Length[A]; 
        Table[Minors[A, i][[1, 1]], {
            i, k}]]; gradf[f_, v_List] := 
      Module[{n, i}, n = Length[v]; 
        Table[(∂_(v[[i]]))f, {
            i, 1, n}]];

arc2d[f_, s_List, r_List, dh___] := 
      Module[{chrds, funcx, arcx}, 
        chrdx = Graphics[{Dashing[{0.02}], Thickness[0.008], 
              Line[{{r[[2]], 
                    f /. r[[1]] -> 
                        r[[2]]}, {
                    r[[3]], 
                    f /. r[[1]] -> 
                        r[[3]]}}]}]; 
        If[dh == Null, chrdx = Graphics[{}]]; 
        funcx = Plot[f, s, DisplayFunction -> Identity]; 
        arcx = Plot[f, r, PlotStyle -> {Thickness[0.009]}, 
            DisplayFunction -> Identity]; {chrdx, funcx, arcx}];

tangentline[f_, s_List, r_] := 
      Block[{tx, $DisplayFunction = Identity}, 
        tx = (f /. s[[1]] -> r)
               + ((∂_(s[[1
                        ]]))f /. 
                  s[[1]] -> r) 
              (s[[1]] - r); {
          Graphics[{}], 
          Plot[{f, tx}, s, PlotStyle -> {{}, {Thickness[0.008]}}]}];

fjohn[a_List, constraints_List, object_, v_List, label___]:=Module[{temp, n, i, j, rul1, m, vlist, glist, vectf, tempf, index, nl, 
              nw},Clear[i, n, rul1, x, y];nl = True;nw = False;If[label == "cone", nw = True];If[label == Null, {nl = False, nw = False}];rul1 = {
                  v[[1]] -> 
                    a[[1]], 
                  v[[2]] -> 
                    a[[2]]};n = Length[constraints];binding = {};glist = {};vlist = {};index = {};Do[If[Chop[
                      N[constraints[[i
                            ]] /. rul1]] == 0, {
                    binding = 
                      Append[binding, 
                        constraints[[i
                          ]]], 
                    index = Append[index, i]}], {i, 1, n}];tempf = gradf[object, v] /. rul1;tempf = (tempf/(Sqrt[tempf . tempf]));If[nw, tempf = -tempf];vectf = 
                vector2d[tempf, 
                  If[nl, If[nw, "-grad f", "grad f"], ""]];m = Length[binding];Do[vlist = 
                  Append[vlist, 
                    gradf[binding[[k]], 
                        v] /. rul1], {k, 1, m}];Do[glist=Append[glist,vector2d[vlist[[j
                                  ]]/(Sqrt[vlist
                                      [[j
                                      ]] . 
                                    vlist[[j
                                      ]]]),If[nl,ToString["grad g\!\(index \\[LeftDoubleBracket] j \\[RightDoubleBracket]\)"],""]]],{j, 1, m}];glist = Append[glist, vectf];glist];

Clear[separableode]; separableode[x_, y_, P_, Q_] := 
      Module[{z, w, Pz, Qw, solz, solw, testz, testw, lst}, testz = True; 
        testw = True; Pz = P /. x -> z; Qw = Q /. y -> w; 
        gensol = 
          ∫Pz \[DifferentialD]z + 
              ∫Qw \[DifferentialD]w == C[1]; 
        solz = Solve[gensol, z]; 
        If[solz == {} || Head[solz] == Solve, testz = False]; 
        solw = Solve[gensol, w]; lst = {gensol}; 
        If[solw == {} || Head[solw] == Solve, testw = False]; 
        If[testz, lst = Append[lst, solz]]; 
        If[testw, lst = Append[lst, solw]]; 
        lst = lst /. {z -> x, w -> y}];

linearode[x_, P_, Q_] := 
      Module[{z, mu, temp, Pz, Qz}, Pz = P /. x -> z; 
        Qz = Q /. x -> z; mu = Exp[∫Pz \[DifferentialD]z]; 
        temp = ((∫mu Qz \[DifferentialD]z + C[1])/mu); 
        temp /. z -> x]; tangentfield[f_, {x_, a_, b_}, {y_, c_, d_}, n_Integer, opt___] := 
      Module[{incx, incy, slope, pt, del, i, j, m, lineset, ptset}, {
          incrx = ((b - a)/(n - 1)); incry = ((d - c)/(n - 1)); 
          del = (incrx/10); slope = {}; pt = {}; 
          Do[Do[pt = Append[pt, {a + i incrx, c + j incry}], {i, 0, 
                n - 1}], {j, 0, n - 1}]; 
          Do[Do[slope = 
                Append[slope, 
                  f /. {x -> a + i incrx, y -> c + j incry}], {i, 
                0, n - 1}], {j, 0, n - 1}]; lineset = {}; ptset = {}; 
          Do[{end1 = {pt[[m, 1]] - del, 
                  pt[[m, 2]] - 
                    slope[[m]] del}, 
              end2 = {pt[[m, 1]] + del, 
                  pt[[m, 2]] + 
                    slope[[m]] del}, 
              lineset = Append[lineset, Line[{end1, end2}]], 
              ptset = Append[ptset, 
                  Point[pt[[m]]]]}, {m, 
              1, (n^2)}]; 
          Graphics[{Thickness[0.0075], lineset, PointSize[0.01], ptset}, 
            opt]}]; \(\(wronskian[lst_List, x_] := 
      Module[{k, i}, 
        Det[Table[∂\_{x, k}lst, {k, 0, Length[lst] - 1}]]]; \)\)\(\(wis[F_, lst_List, x_] := 
      Module[{n, k, i, col, w, an, wdet}, an = 1; n = Length[lst]; wdet = {}; 
        w = Table[∂\_{x, kk}lst, {kk, 0, Length[lst] - 1}]; 
        col = Table[If[i ≠ n, 0, F\/an], {i, 1, n}]; 
        Do[wdet = 
            Append[wdet, Det[Transpose[ReplacePart[Transpose[w], col, k]]]], {
            k, 1, n}]; wdet]; \)\)

Clear[trajectory, yy, ran]; trajectory[yy_List, ran_List, opts___] := 
      Module[{u, len, t, ptlst, maxdis, maxpos, n1, plt}, 
        plt = ParametricPlot[{yy[[1]], 
              yy[[2]]}, ran, 
            DisplayFunction -> Identity, opts]; 
        u = Position[plt, Line][[1
            ]]; len = Length[u]; 
        u[[len]] = 
          u[[len]] + 1; 
        t = Prepend[u, plt]; ptlst = Part@@t; len = Length[ptlst]; 
        maxdis = 0; maxpos = 1; 
        For[i = 1, i < len, i++, 
          If[(ptlst[[i]] - 
                    ptlst[[i + 1]]) . 
                (ptlst[[i]] - 
                    ptlst[[i + 1]]) > 
              (maxdis^2), maxpos = i]]; 
        If[maxpos == len - 1, maxpos = maxpos - 1]; 
        If[maxpos == 1, maxpos = maxpos + 2]; n1 = maxpos - 1; 
        Show[plt, 
          Graphics[
            Arrow[ptlst[[n1]], 
              ptlst[[n1 + 1]]]], 
          DisplayFunction -> $DisplayFunction]; ];

End[]; Protect[fibprime, directed, vector2d, directed3d, vector3d, 
      addition2d, subtraction2d, projection2d, projection3d, rotation2d, 
      rank, rowop, colop, submatrix, minor, cofactor, aug, colrep, singQ, 
      signQL, borderB, gradf, hessian, conhess, arc2d, tangentline, fjohn, 
      separableode, linearode, tangentfield, wronskian, wis, trajectory]; 
    EndPackage[];