Example 1: The neoclassical model with money in utility function that is separable is money and consumption. The economy experiences a sustained rise in productivity. One solves for the global nonlinear saddle path in two stages. In the first stage, one solves the (i,k,z) system by reverse shooting in two dimesion in (k,z)-space. Reverse shooting is possible here as this subsystem does not have a unit root. In the second stage, one solves the (c,b) system with unit root by forward shooting. In this model, both c and m jump right away to their final steady state values.
THE OUTPUT OF THIS FILE CORRESPONDS TO FIGURE 2 OF THE PAPER.

Off[General::spell]; Off[General::spell1]; Off[Unset::norep]; 
    Off[FindMinimum::lstol]; Off[FindMinimum::reged]; 
    Off[InterpolatingFunction::dmval]; Off[Precision::precsm];

(* 
          This cell sets the precision of calculations. *) 
{precision  =  
            100, $MinPrecision = digits = precision};

Initialization

SetDirectory["c:"];

{τ = (1/4), ρ = (1/10), δ = (1/20), σ = 
            (1/2), Ω = 3}//N

{(Subscript[bbyY, o]) = -(25/100), 
        (Subscript[mbyY, o]) = (15/200), (Subscript[π, o]) = (1/2), (Subscript[θ, ko]) = (1/2)}//N

{(Subscript[z, o]) = 1, (Subscript[c, o]) = 1}//N

{(Subscript[z, f]) = (6/5) (Subscript[z, o]), s = (1/5)}//N

Steady State

U[c_] = (1/(1 - (1/τ))) 
            (c^(1 - (1/τ))) ;(Subscript[U, 1])[c_] = 
          D[U[c], c] ;ϕ[
            m_] = ((Subscript[a, 0])/(1 - 
                  (1/τ))) (((m))^(1 - (1/τ))) ;(Subscript[ϕ, 1])[m_] = 
          D[ϕ[m], m] ;
F[k_, z_] = 
          z ((((Subscript[a, 1]) + (Subscript[a, 2]) (k^(((σ - 1)/σ)))))^((σ/(σ - 1)))) ;(Subscript[F, 1])[k_, z_] = 
          D[F[k, z], k] ;ψ[
            x_] = (v/2) (((x - δ))^2) ;v = 
          (1/(Ω δ)) ;

(Subscript[m, o]) =  (Subscript[mbyY, o]) (Subscript[Y, o]); 
      (Subscript[b, o]) = (Subscript[bbyY, o]) (Subscript[Y, o]); 
      g = (Subscript[π, o]) 
          (Subscript[m, o]); 
eq1 := (Subscript[Y, o]) (1 + ρ (Subscript[bbyY, o])) == (Subscript[c, o]) + δ (Subscript[k, o]) ; 
      eq2 := (((ρ + δ) (Subscript[k, o]))/(Subscript[Y, o])) == (Subscript[θ, ko]); {(Subscript[Y, o]), (Subscript[k, o])} =. ; {(Subscript[Y, o]), (Subscript[k, o])} = {(Subscript[Y, o]), (Subscript[k, o])} /. 
          FindRoot[{ eq1, eq2}, {(Subscript[Y, o]), 1}, {(Subscript[k, o]), 1}, 
            MaxIterations -> 1000, AccuracyGoal -> $MinPrecision, 
            WorkingPrecision -> 
              digits]; 
(Subscript[a, 0]) =. ;

      {(Subscript[a, 0])} = {(Subscript[a, 0])}  /.  
          FindRoot[{ (Subscript[ϕ, 1])[(Subscript[m, o])] ==  
                (Subscript[U, 1])[(Subscript[c, o])] (ρ + (Subscript[π, o]))}, {(Subscript[a, 0]), (10^(-3))}, 
            MaxIterations -> 1000, AccuracyGoal -> $MinPrecision, 
            WorkingPrecision -> digits]; {(Subscript[a, 1]), (Subscript[a, 2])} =. ; {(Subscript[a, 1]), 
          (Subscript[a, 2])} = {(Subscript[a, 1]), (Subscript[a, 2])} /. 
          FindRoot[{(Subscript[F, 1])[(Subscript[k, o]), 
                  (Subscript[z, o])]  ==  (ρ + δ), (Subscript[θ, ko]) == (((ρ + δ) (Subscript[k, o]))/F)[(Subscript[k, o]), (Subscript[z, o])]}, {(Subscript[a, 1]), 1}, {(Subscript[a, 2]), 
              1}, MaxIterations -> 1000, AccuracyGoal -> digits, 
            WorkingPrecision -> digits]; 
{{(Subscript[m, o]), 
              (Subscript[c, o]), (Subscript[k, o]), (Subscript[Y, o]), (Subscript[b, o]), g, (Subscript[z, o]), (Subscript[π, o]), (Subscript[c, o])/(Subscript[m, o]), 
              (Subscript[m, o])/(Subscript[Y, o])}} // N // MatrixForm

      {(Subscript[a, 0]), (Subscript[a, 1]), (Subscript[a, 2])} // N

(Subscript[k, f]) =. ; 
      (Subscript[k, f]) = (Subscript[k, f]) /. 
          FindRoot[{(Subscript[F, 1])[(Subscript[k, f]), (Subscript[z, f])] == ρ + δ}, {(Subscript[k, f]), (Subscript[k, o])}, 
            MaxIterations -> 1000, AccuracyGoal -> digits, 
            WorkingPrecision -> digits]; 
      (Subscript[Y, f]) = F[(Subscript[k, f]), 
          (Subscript[z, f])]; 
{{(Subscript[m, o]), (Subscript[c, o]), (Subscript[k, o]), (Subscript[Y, o]), (Subscript[b, o]), g, 
              (Subscript[z, o]), (Subscript[π, o]), (Subscript[c, o])/(Subscript[m, o]), (Subscript[m, o])/(Subscript[Y, o])}, {(Subscript[m, f]), (Subscript[c, f]), (Subscript[k, f]), (Subscript[Y, f]), 
              (Subscript[b, f]), g, (Subscript[z, f]), (Subscript[π, f]), (Subscript[c, f])/(Subscript[m, f]), (Subscript[m, f])/(Subscript[Y, f])}} // N // 
        MatrixForm

Solving for the Global Nonlinear Saddlepath

SOLVING THE PRODUCTION SIDE BY REVERSE SHOOTING (ONE DIMENSIONAL SEARCH ON ϵ-RING IN (k,z)-space)

idot = (k/(ψ''))[i/k]((ρ + δ) (1 + ψ'[
                          i/k])+ψ[i/k]-(Subscript[F, 1])[k, z]);kdot = i - δ k;zdot = s ((Subscript[z, f]) - z);sysvarsubs={i->i[t],k->k[t],z->z[t]};system={idot,kdot,zdot};

sys=Thread[{i'[t],k'[t],z'[t]}==(-system /. sysvarsubs),List];ini:=Thread[{i[0],k[0],z[0]}=={iin,kin,zin},List];
  
box:=  i[t]<((Min[δ (Subscript[k, o]),δ (Subscript[k, f])])/10)||i[t]>2Max[δ (Subscript[k, o]),δ (Subscript[k, f])]|| k[t]<((Min[(Subscript[k, o]),(Subscript[k, f])])/10)||k[t]>2Max[(Subscript[k, o]),(Subscript[k, f])]|| z[t]<((Min[(Subscript[z, o]),(Subscript[z, f])])/(                         10))||z[t]>2Max[(Subscript[z, o]),(Subscript[z, f])];
  
printandplot:=(Print["Extreme Left Trajectory = ", 
                  j1, "   Closest Trajectory = ", 
                  jmin, "   Extreme Right Trajectory = ", 
                  j2, "    No. of Subdivisions of the ϵ-ring = \
", ntraj];Print["Current Minimum Distance = ", 
                  gmindis, "  Time of Closest Approach = ", gtmin];
ParametricPlot[{{k[t] ,z[t]}/.gminsol,{(Subscript[k, o])+(((Subscript[k, f])-(Subscript[k, o]))/100)Sin[
                                  t],(Subscript[z, o])+(((Subscript[z, f])-(Subscript[z, o]))/100) Cos[t]}},{t, 0, gtmin},PlotRange -> All];)  ;
  
sysdiff[t1_, t2_]:={Max[Abs[Table[((i[gtmin - t]/.gminsol)-(i[gtmin1 - t]/.gminsol1))[[1]],{t, t1, t2, 1/10}]]],Max[Abs[Table[((k[gtmin - t]/.gminsol)-(k[gtmin1 - t]/.gminsol1))[[1]],{t, t1, t2, 1/10}]]],Max[Abs[Table[((z[ gtmin - t]/.gminsol)-(z[gtmin1 - t]/.gminsol1))[[1]],{t, t1, t2, 1/10}]]]};

distance[sol_, tsol_]:=If[firstiter == 1,
Block[{dis, dislast, tdis, steps, tbeg, tend, intsz},dis =  (10^7);dislast = (10^8);steps = Max[20, 4 Floor[tsol]] ;tbeg = 0;tend = tsol;intsz = tend - tbeg;
While[ (dislast - dis) > tol (10^(-5)),dislast = dis;traj=Table[Norm[Flatten[{k[t ]-(Subscript[k, o]),z[t]-(Subscript[z, o ])}/.sol],2],{t, tbeg, tend, (intsz/steps)}];
dis = Min[traj];tdis = 
                            tbeg + (Position[traj, dis][[1, 1]] - 
                                    1) (intsz/steps);tbeg = Max[tdis - (intsz/steps), 0];tend = Min[tdis + (intsz/steps), tsol];intsz = tend - tbeg;steps = 20];{dis, tdis}],
Block[{dis, tdis},tbeg = 0;tend = tsol;min=FindMinimum[First[Norm[Flatten[{k[t]-(Subscript[k, o]),z[t]-(Subscript[z, o])} ]]/.sol] ,{t, 9  tsol/10, tbeg, tsol},WorkingPrecision -> digits,AccuracyGoal -> 5 - Log[10, tol]];{dis, tdis} = {min[[1]], 
                        min[[2, 1, 2]]}]];

(* MESSAGES - MESSAGES. *)
diffprec := 
            Print["Difference in values of variables due to latest increase \
in precision = ",  sysdiff[0, T] // N];
diffeps := 
            Print["Difference in values of variables due to latest decrease \
in epsilon = ", sysdiff[0, T] // N];
satisprec= "THE PRECISION IS SATISFACTORY. The variables changed by less\! than\! 'tol' due to latest increase in 'precision' for current ϵ-ring.";
notsatisprec="The precision is NOT satisfactory. The variables changed by greater than\! 'tol' due to latest increase in 'precision' for current ϵ-ring.";
notsatiseps="The epsilon rings are NOT tight enough. The variables changed by greater than\! 'tol' due to latest decrease in 'epsilon'.";
increaseprec:=Print["Increased 'precision' by 4 digits to ",precision," digits."];decreaseps:=Print["Decreased 'epsilon' by a half to ",N[epsilon, precision],"."];
step2="Starting Step 2: Checking for the precision of the candidate solution.";step3="Starting Step 3: Checking for the accuracy of the candidate solution.";
congratulations= "CONGRATULATIONS! The variables changed by less\! than\! 'tol' due to latest decrease in 'epsilon'.";
done = "You are done with solving the PRODUCTION sub-system! Continue with the next cell.";

(* INITIALIZATION CELL - Initialization for finding the required trajectory. *)
tol=(10^(-5));                                                  ;
T = 20;                 
epsilon= (10^(-3))Norm[{(Subscript[k, f])-(Subscript[k, o]),(Subscript[z, f])-(Subscript[z, o])}];      ;
precisionini=16;                             ;

(* MAIN CELL FOR SOLVING THE PRODUCTION SUB - SYSTEM *)
Timing[initialize;Tmax = 4000;firstring  =  01;firstrun  =  01;precision = precisionini;
While[True,  
Label[precbegin];$MinPrecision = digits = precision;ntraj = 20;j1 = 0;j2 = ntraj;gmindis = (10^100);firstiter = 1;Label[begin];
For[j = j1,j ≤  j2,j++,theta = 2  π j/ntraj;{dkin=epsilon Cos[theta],dzin=epsilon Sin[theta]};{iin,kin,zin}={δ (Subscript[k, f]),(Subscript[k, f])+dkin,(Subscript[z, f])+dzin};
If[((IntegerPart[j/2] - j/2) ≠ 
                                0 || firstiter == 1),If[firstiter == 0&&Max[Abs[{kin/pmin[[1]]-1,zin/pmin[[2]]-1}]]<(10^(-precision)),precision = precision + 4;Print["Suspecting loss in precision."];Print[increaseprec];Goto[precbegin]];
sol=NDSolve[SetPrecision[{sys, ini}, 
                                    precision],{i, k, z},{t, 0, Tmax},\(Method -> {"EventLocator", \ "Event" -> box, "EventAction" :> \ 
                                        Throw[
                                        Null, \ "StopIntegration"], "EventLocationMethod" -> "StepBegin"}\),MaxSteps -> 100000,WorkingPrecision -> digits];tsol = sol[[1, 1, 2, 1, 1, 2]];{dis, tdis} = distance[sol, tsol];
If[firstiter == 1, 
                              Print["Trajectory = ", 
                                j, "     Duration of Solution = ", 
                                tsol // N, "     Minimum Distance = ", 
                                dis, "     Time of Closest Approach = ", 
                                tdis]];
If[dis ≤ gmindis,gmindis = dis;gtmin = tdis;jmin = j;gminsol = sol ;pmin={kin,zin};If[gmindis < tol/3,Print["Target Tolerance = \
", tol/3 // N, "   Tolerance Achieved = ", gmindis]; Goto[end] ];];];];firstiter = 0;printandplot;{j1, j2} = 2 {jmin - 1, jmin + 1};{ntraj, jmin} = 2 {ntraj, jmin};Goto[begin];Label[end];
  
Which[firstring  ==  1  &&  
                      firstrun  ==  1,If[gtmin < T + 60, epsilon = epsilon/2; 
                      decreaseps, {gminsol, gtmin} >> solution1; 
                      Print[gtmin]; Print[step2]; precision = precision + 4; 
                      increaseprec; firstrun = 0],firstring  ==  1  &&  
                      firstrun  ==  0,{gminsol1, gtmin1} = << solution1; {gminsol, 
                        gtmin} >> solution1; diffprec; 
                    If[Max[sysdiff[0, T]] > tol, Print[notsatisprec]; 
                      precision = precision + 4; increaseprec, 
                      Print[satisprec]; Print[step3]; firstring = 0; 
                      firstrun = 1; epsilon = epsilon/2; decreaseps],firstring  ==  0   &&  
                      firstrun  ==  1,{gminsol, gtmin} >> solution2;{gminsol1, gtmin1} = << solution1;If[Max[sysdiff[0, T]] > tol,firstrun = 0; precision = precision + 4; 
                          increaseprec,diffeps;Print[congratulations];Print[done];Goto[end1] ],
firstring  ==  
                        0   &&  firstrun  ==  0,{gminsol1, gtmin1} = << solution2;{gminsol, gtmin} >> solution2;diffprec;
If[Max[sysdiff[0, T]] > tol,Print[notsatisprec]; 
                          precision = precision + 4; increaseprec,Print[satisprec];{gminsol1, gtmin1} = << solution1;diffeps;If[Max[sysdiff[0, T]] > tol,{gminsol, gtmin} >> solution1; 
                                Print[notsatiseps]; firstrun = 1; 
                                epsilon = epsilon/2; decreaseps,Print[congratulations];Print[done];Goto[end1]]]];];Label[end1];
  
IKprec = Precision[gminsol];For[j = T + 1,j ≤ IntegerPart[gtmin1],j++,If[Max[sysdiff[j - 1, j]] > tol,Print["The solution is accurate within tolerance 'tol' upto ",j - 1," years. Continue with next cell."];Break[]]];If[j == IntegerPart[gtmin1] + 1,Print["The solution is accurate upto ",j - 1," years. Continue with next cell."]]]

SOLVING THE CONSUMPTION-BOND SYSTEM BY FORWARD SHOOTING (ONE DIMENSIONAL SEARCH IN c-space)

bdot[t_?NumberQ] := First[
                F[k[gtmin - t], 
                      z[gtmin - t]]  - (i[
                          gtmin - t] + ψ[i[gtmin - t]/k[gtmin - t]] 
                          k[gtmin - t]) /. gminsol];
  
sysdiff[t1_, t2_] := 
            Max[Abs[Table[(b[t] /. bminsol) - (b[t] /. bminsol1), {t, 
                    t1, t2, 1/10}]]];

Clear[distance];
distance[sol_, tsol_]:=
Block[{dis, dislast, tdis, steps, tbeg, tend, intsz},dis =  (10^7);dislast = (10^8);steps = Max[20, 4 Floor[tsol]] ;tbeg = 5/10  tsol;tend = tsol;intsz = tend - tbeg;
While[ (dislast - dis) > 
                        Max[tol, bmindis] (10^(-5)),dislast = dis;traj=Table[Norm[Flatten[{b[t]-(Subscript[b, f])}/.sol],2],{t, tbeg, tend, (intsz/steps)}];
dis = Min[traj];tdis = 
                          tbeg + (Position[traj, dis][[1, 1]] - 
                                  1) (intsz/steps);tbeg = Max[tdis - (intsz/steps), 0];tend = Min[tdis + (intsz/steps), tsol];intsz = tend - tbeg;steps = 20];{dis, tdis}];

(* MESSAGES - MESSAGES. *)
\(difftmin := 
            Print[\*""Difference in values of variables between two paths with (Subscript[t, min])s differing by 5 years = "", sysdiff[0, T]];\)
notsatistmin="The (Subscript[(t), min]) is NOT large enough. The variables changed by greater than\! 'tol' due to latest increase in (Subscript[(t), min]).";
satistmin= "THE ACCURACY CHECK PASSED SUCCESSFULLY. The variables changed by less\! than\! 'tol' due to latest increase in (Subscript[(t), min]).";
diffprec := 
            Print["Difference in values of variables due to latest increase \
in precision = ",  sysdiff[0, T]];
notsatisprec="The precision is NOT satisfactory. The variables changed by greater than\! 'tol' due to latest increase in 'precision'.";
satisprec="THE PRECISION IS SATISFACTORY AS WELL. The variables changed by less\! than\! 'tol' due to latest increase in 'precision'.";increaseprec:=Print["Increased 'precision' by 4 digits to ",precision," digits."];
manualincprec:=Print["Increase precision by 4 in the INITIALIZATION CELL, quit Mathematica kernal, and re-execute the file from beginning."];
step2="Starting Step 2: Checking for the accuracy of the candidate solution.";step3="Starting Step 3: Checking for the precision of the candidate solution.";
congratulations= "CONGRATULATIONS! You are done with solving CONSUMPTION-BOND UNIT ROOT system!";

(* 
          Set up for solving for the global nonlinear saddle path *) 
Clear[compdis];
compdis[p_, ϵ_]:=Block[{bpmin, epsilon},iflag = 0;it = 1;{bpmin, epsilon} = {p, ϵ};Label[begin];Print["\nBEGINNING ITERATION ",it," WITH CENTER OF ϵ-segment at ",SetPrecision[bpmin, precision]];
For[j = -1,j ≤  1  ,j = j + 2,{cin,bin}={bpmin + j epsilon,(Subscript[b, o])};(Subscript[c, f]) = cin;(Subscript[b, f]) = ((Subscript[c, f]) + δ (Subscript[k, f]) - 
                                  (Subscript[Y, f]))/ρ;sol = First[NDSolve[SetPrecision[{b'[t] == 
                                        bdot[t] + ρ b[t] - cin, 
                                        b[0] ==  bin},precision],{b},{t, 0, Tmax},WorkingPrecision -> digits ,MaxSteps -> 100000]];{dis, tdis} = 
                          distance[sol, sol[[1, 2, 1, 1, 2]]];jdis[j] = dis;jtdis[j] = tdis;dissol[j] = sol;jp[j]=cin;jb[j]=b[sol[[1, 2, 1, 1, 2]]]/.sol;jbf[j]=(Subscript[b, f]);];If[
                    jbf[1] < jb[1]  &&  jbf[-1] > jb[-1], 
                    Print["ERROR"]; manualincprec; Abort[]];
If[
                    jbf[1] < jb[1]  &&  
                      jbf[-1] < jb[-1], {bmindis, btmin, bminsol, 
                        bpmin} = {jdis[1], jtdis[1], dissol[1], jp[1]}, 
                    If[jbf[1] > jb[1]  &&  
                        jbf[-1] > jb[-1], {bmindis, btmin, bminsol, 
                          bpmin} = {jdis[-1], jtdis[-1], 
                          dissol[-1], jp[-1]}, 
                      epsilon = epsilon/2]];
Print["AT THE END OF ITERATION = ",it,"     TIME OF CLOSEST APPROACH = ",btmin];
If[firstsim == 1,If[btmin > T + 5  &&  iflag == 0, 
                          iflag = 
                            1; {bminsol, btmin, bmindis, bpmin, 5  epsilon} >> 
                            bsolution1; Print[""]; Print[step2];];
If[btmin > 
                                Max[T + 5 , (<< bsolution1)[[2]] + 
                                    5] &&  iflag == 1 ,{bminsol1, 
                                  btmin1} = (<< bsolution1)[[{1, 
                                      2}]];If[sysdiff[0, T] > tol,Print[""]; difftmin; 
                                  Print[
                                    notsatistmin]; {bminsol, btmin, bmindis, 
                                      bpmin, 5 epsilon} >> bsolution1,Print[""];difftmin;Print[satistmin];{bminsol, btmin, bmindis, bpmin, 
                                        5 epsilon} >> bsolution1;Print[step3];Goto[end];]],
If[btmin > (<< bsolution2)[[2]] &&  
                              bmindis < (<< bsolution2)[[3]],{bminsol, btmin, bmindis, bpmin, 
                                  5  epsilon} >> bsolution1;Goto[end];];];
If[(epsilon/bpmin) < (10^(-precision)),Print["You may be losing precision in calculations."];gtmin = (10^100);Goto[end];]it++;Goto[begin];Label[end];{bminsol, btmin, bmindis, 
                    bpmin, 5 epsilon}];

Timing[precision = precisionini;firstsim = 1;Tmax = gtmin;bmindis = 1;
While[True,Label[begin];$MinPrecision = digits = precision;If[firstsim == 0 ,{pini, 
                      epsilon} = (<< bsolution2)[[{4, 5}]],{pini, epsilon}={(Subscript[c, o]), (Subscript[c, o])/4}];
  ;
{bminsol, btmin, bmindis, bpmin, eps} = 
                compdis[pini, epsilon];If[gtmin == (10^100),precision = precision + 4;;Goto[begin]];
   
If[firstsim  ==  1,firstsim = 0; 
                  If[precision ==  IKprec, Print["ERROR"]; 
                    manualincprec; Abort[]]; precision = precision + 4; 
                  increaseprec; {bminsol, btmin, bmindis, bpmin, eps} >> 
                    bsolution2,{bminsol1, btmin1, 
                        bpmin1} = (<< bsolution2)[[{1, 2, 4}]];If[Max[sysdiff[0, T]] > tol,Print[""]; diffprec; Print[notsatisprec]; 
                        If[precision == IKprec, Print["ERROR"]; 
                          manualincprec; Abort[]]; precision = precision + 4; 
                        increaseprec; {bminsol, btmin, bmindis, bpmin, eps} >> 
                          bsolution2,Print[""];diffprec;Print[satisprec];Print[""];Print[congratulations];Goto[end]]];];
  ;
Label[end];For[j = T + 1,j ≤ IntegerPart[Min[btmin1, btmin]],j++,
If[Max[sysdiff[j - 1, j]] > tol,Print["The solution is accurate within tolerance 'tol' upto ",j - 1," years. Continue with next cell."];Break[]]];If[j == IntegerPart[Min[btmin1, btmin]] + 1,Print["The solution is accurate upto ",j - 1," years. Continue with next cell."]]; ]

RESULTS

{gminsol, 
            gtmin} = (<< solution2)[[{1, 2}]]; {bminsol, 
            bpmin} = (<< bsolution1)[[{1, 4}]];
(Subscript[c, f]) = bpmin; 
          (Subscript[b, f]) = (((Subscript[c, f]) + δ (Subscript[k, f]) - (Subscript[Y, f]))/ρ); {(Subscript[m, f]), (Subscript[π, f])} =. ;
{(Subscript[m, f]), (Subscript[π, f])} = {(Subscript[m, f]), (Subscript[π, f])}  /.  
                FindRoot[{ (Subscript[ϕ, 1])[(Subscript[m, f])] ==  
                      (Subscript[U, 1])[(Subscript[c, f])] (ρ + (Subscript[π, f])), 
                    g == (Subscript[π, f]) (Subscript[m, f])}, {(Subscript[m, f]), 
                    (Subscript[m, o])}, {(Subscript[π, f]), (Subscript[π, o])}, MaxIterations -> 1000, 
                  AccuracyGoal -> $MinPrecision, 
                  WorkingPrecision -> digits];
{{(Subscript[m, o]), (Subscript[c, o]), (Subscript[k, o]), (Subscript[Y, o]), (Subscript[b, o]), g, (Subscript[z, o]), (Subscript[π, o]), 
                  (Subscript[c, o])/(Subscript[m, o]), (Subscript[m, o])/(Subscript[Y, o])}, {(Subscript[m, f]), (Subscript[c, f]), (Subscript[k, f]), (Subscript[Y, f]), (Subscript[b, f]), g, 
                  (Subscript[z, f]), (Subscript[π, f]), (Subscript[c, f])/(Subscript[m, f]), (Subscript[m, f])/(Subscript[Y, f])}} // N // 
            MatrixForm

tplot =T;
pz=Plot[{If[t < 0,(Subscript[z, o]),z[gtmin - t]/.gminsol]},{t, -2, tplot},PlotRange -> { .95, 1.21},AxesOrigin -> {-2,  .95},AxesLabel -> {"t", z},PlotStyle -> {{Thickness[ .007], 
                  Dashing[{ .02, 0}]}},DisplayFunction -> Identity];pc=Plot[{If[t < 0,(Subscript[c, o]),(Subscript[c, f])]},{t, -2, tplot},PlotRange -> { .95, 1.21},AxesOrigin -> {-2,  .95},AxesLabel -> {"t", c},PlotStyle -> {{Thickness[ .007], 
                  Dashing[{ .02, 0}]}},DisplayFunction -> Identity];
pm=Plot[{If[t < 0,(Subscript[m, o]),(Subscript[m, f])]},{t, -2, tplot},PlotRange -> { .088,  .115},AxesOrigin -> {-2,  .088},AxesLabel -> {"t", m},PlotStyle -> {{Thickness[ .007], 
                  Dashing[{ .02, 0}]}},DisplayFunction -> Identity];pπ=Plot[{If[t < 0,(Subscript[π, o]),(Subscript[π, f])]},{t, -2, tplot},PlotRange -> { .40,  .5},AxesOrigin -> {-2,  .4},AxesLabel -> {"t", π},PlotStyle -> {{Thickness[ .007], 
                  Dashing[{ .02, 0}]}},DisplayFunction -> Identity];
pi=Plot[{If[t < 0,δ (Subscript[k, o]),i[gtmin - t]/.gminsol]},{t, -2, tplot},PlotRange -> { .2,  .265},AxesOrigin -> {-2,  .2},AxesLabel -> {"t", "I"},PlotStyle -> {{Thickness[ .007], 
                  Dashing[{ .02, 0}]}},DisplayFunction -> Identity];pk=Plot[{If[t < 0,(Subscript[k, o]),k[gtmin - t]/.gminsol]},{t, -2, tplot},PlotRange -> {4.05, 4.81},AxesOrigin -> {-2, 4.05},AxesLabel -> {"t", K},PlotStyle -> {{Thickness[ .007], 
                  Dashing[{ .02, 0}]}},DisplayFunction -> Identity];
pb=Plot[{If[t < 0,(Subscript[b, o]),b[t]/.bminsol]},{t, -2, tplot},PlotRange -> {-2, 0},AxesOrigin->{-2,Automatic},AxesLabel -> {"t", b},PlotStyle -> {{Thickness[ .007], 
                  Dashing[{ .02, 0}]}},DisplayFunction -> Identity];pY=Plot[{If[t < 0,(Subscript[Y, o]),F[k[gtmin - t], z[gtmin - t]]/.gminsol]},{t, -2, tplot},PlotRange -> {1.2, 1.6},AxesOrigin -> {-2, 1.2},AxesLabel -> {"t", Y},PlotStyle -> {{Thickness[ .007], 
                  Dashing[{ .02, 0}]}},DisplayFunction -> Identity];
Show[GraphicsArray[{{pz,pc},{pm,pπ},{pi,pk},{pY,pb}}]]

tplotNL = 
        6;
color = {{RGBColor[1, 0, 
              0]}, {RGBColor[0, 1, 0]}};
pCxall=Plot[{If[t < 0,0,((Cx[(Subscript[λ, 1, f]),DD[t],t]-(Subscript[C, o]))/((Subscript[C, o])))100/.gminsol]},{t, -1, tplotNL},PlotRange -> All,AxesLabel -> {"t", ""},PlotLabel -> "Non-durables Consumption",AxesOrigin -> {-1, 0},PlotStyle -> {{Thickness[ .007], 
                  Dashing[{ .05, 0}]}},
DisplayFunction -> Identity];
pSxall=Plot[{If[t < 0,0,((Sx[(Subscript[λ, 1, f]),λ2[t],DD[t],t]-(Subscript[S, o]))/((Subscript[S, o])))100/.gminsol]},{t, -1, tplotNL},PlotRange -> All,AxesLabel -> {"t", ""},PlotLabel -> "Durables Spending",AxesOrigin -> {-1, 0},PlotStyle -> {Thickness[ .007], Dashing[{ .05, 0}]},DisplayFunction -> Identity];
pExall=Plot[{If[t < 0,0,((Sx[(Subscript[λ, 1, 
                                    f]),λ2[t],DD[t],t]+Cx[(Subscript[λ, 1, 
                                    f]),DD[t],t]-(Subscript[S, o])-(Subscript[C, o]))/((Subscript[S, o]) + (Subscript[C, o])))100/.gminsol],If[t < 0 , 0, 
                  If[t < 3, 
                    0.03441110832850172` 100, -0.00893013093579049` 100]]},{t, -1, tplotNL},PlotRange -> All,AxesLabel -> {"t", ""},PlotLabel -> "Total Consumption",AxesOrigin -> {-1, 0},PlotStyle -> {{Thickness[ .007], 
                  Dashing[{ .05, 0}]}, {Thickness[ .007], 
                  Dashing[{ .02,  .02}]}},DisplayFunction -> Identity];
pZall = Plot[{If[t < 0,0,((Z[t] - (Subscript[Z, o]))/(Subscript[Q, o]))/.gminsol]},{t, -1, tplotNL},PlotRange -> All,AxesLabel -> {"Years", ""},AxesOrigin -> {-1, 0},PlotLabel -> "Total Net Foreign Assets",PlotStyle -> {Dashing[{ .05, 0}], 
                Dashing[{ .02,  .02}]},DisplayFunction -> Identity];
Show[
      GraphicsArray[{pExall,   pCxall, 
          pSxall}]  ]

Labels={"m","c","i","b","k","z","Y","pi"};InitialSteadyState={(Subscript[m, o]),(Subscript[c, o]),δ (Subscript[k, o]),(Subscript[b, o]),(Subscript[k, o]),(Subscript[z, o]),(Subscript[Y, o]),(Subscript[π, o])}//N;FinalSteadyState={(Subscript[m, f]),(Subscript[c, f]),δ (Subscript[k, f]),(Subscript[b, f]),(Subscript[k, f]),(Subscript[z, f]),(Subscript[Y, f]),(Subscript[π, f])}//N;
(Dat=Table[Flatten[{t,(Subscript[m, f]),(Subscript[c, f]),i[gtmin - t],b[t],k[gtmin - t],z[gtmin - t],F[k[gtmin - t], z[gtmin - t]],(Subscript[π, f])}/.bminsol/.gminsol],{t,0,T,1/4}])                       //MatrixForm;
(Dat = 
                Join[{Join[{"Initial Steady State"}, InitialSteadyState], 
                    Join[{"Final Steady State"}, FinalSteadyState], 
                    Join[{Time}, Labels]}, Dat]) // N // 
          MatrixForm