Example 4: The neoclassical model with money in utility function that is separable is money and consumption. The economy has two sectors with labor completely mobile across sectors. The capital stock in nontraded sector is fixed as in Example 2. The economy experiences a temporary rise in productivity from t = 0 to t = ts. As in Example 2 with permanent shock, in this case as well, one has to solve the production and the consumption blocks together. Finally, one needs to solve for m by forward shooting. However, now one has to take additional care: 
THIS PROGRAM USING THE PENALTY FUNCTION APPROACH. IF ONE IS UNABLE TO SOLVE THE SYSTEM FOR DEFAULT INITIAL VALUES OR SUITABLE GUESS VALUES OF pini THEN ONE WILL NEED TO SOLVE FIRST THE PERMANENT CASE JUST BY SETTING ts TO BE A LARGE NUMBER AND THEN USING THE pmin OBTAINED FROM THAT CASE AS pini FOR THE TEMPORARY CASE. 

THE OUTPUT OF THIS FILE CORRESPONDS TO FIGURE 11 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), α = (1/2), β = (1/2), Ω = 3}//N

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

{(Subscript[z, o]) = 1, (Subscript[Ex, o]) = 1, (Subscript[Pn, o]) = 1, lbar = 1}//N

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

Steady State

c[
          Pn_] = (((((((Subscript[k, 0])))^β) + ((((Subscript[k, 1])))^β) (Pn^(1 - β))))^((1/(1 - β)))); 
      V[Pn_, Ex_] = (1/(1 - 
                (1/τ))) ((((Ex/c)[Pn]))^(1 - (1/τ))); 
      (Subscript[V, 1])[Pn_, Ex_] = D[V[Pn, Ex], Pn]; 
      (Subscript[V, 2])[Pn_, Ex_] = 
        D[V[Pn, Ex], Ex]; 
ϕ[Pn_, 
          m_] = ((Subscript[a, 0])/(1 - 
                (1/τ))) ((((m/c)[Pn]))^(1 - (1/τ))); (Subscript[ϕ, 2])[
          Pn_, m_] = 
        D[ϕ[Pn, m], m]; ψ[x_] = (v/2) (((x - δ))^2); 
      v = (1/(Ω δ)); 
F[kx_, lx_, 
          z_] = z ((((Subscript[a, 1]) (lx^(((σ - 1)/σ))) + (Subscript[a, 2]) (kx^(((σ - 1)/σ)))))^((σ/(σ - 1)))); 
      (Subscript[F, 1])[kx_, lx_, z_] = D[F[kx, lx, z], kx]; 
      (Subscript[F, 2])[kx_, lx_, z_] = 
        D[F[kx, lx, z], lx]; 
G[kn_, ln_, z_] = 
          z ((((Subscript[a, 3]) (ln^(((σ - 1)/σ))) + (Subscript[a, 4]) (kn^(((σ - 1)/σ)))))^((σ/(σ - 1)))) ;(Subscript[G, 1])[kn_, ln_, z_] = 
          D[G[kn, ln, z], kn] ;(Subscript[G, 2])[kn_, ln_, z_] = 
          D[G[kn, ln, z], ln] ;

(Subscript[Pk, o]) = ((((Subscript[Pn, o])))^α); 
      (Subscript[Cn, o]) = (((Subscript[γ, o]) (Subscript[Ex, o]))/(Subscript[Pn, o])); 
      (Subscript[Y, o]) = (Subscript[Pn, o]) (Subscript[Qn, o]) + (Subscript[Qx, o]); (Subscript[ln, o]) = lbar - (Subscript[lx, o]); 
      (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[Ex, o]) + (Subscript[Pk, o]) δ ((Subscript[kx, o]) + 0  (Subscript[kn, o])) ; 
      eq2 :=  (Subscript[Qn, o]) == 
          (Subscript[Cn, o]) + ((α  (Subscript[Pk, o]) δ ((Subscript[kx, o]) + 0  (Subscript[kn, o])))/(Subscript[Pn, o])); 
      eq3 := (((ρ + δ) (Subscript[Pk, o]) (Subscript[kx, o]))/(Subscript[Qx, o])) == (Subscript[θ, kxo]); 
      eq4 := (((ρ + δ) (Subscript[Pk, o]) (Subscript[kn, o]))/((Subscript[Pn, o]) (Subscript[Qn, o]))) == (Subscript[θ, kno]); {(Subscript[Qx, o]), (Subscript[Qn, o]), (Subscript[kx, o]), 
          (Subscript[kn, o])} =. ; 
{(Subscript[Qx, o]), (Subscript[Qn, o]), (Subscript[kx, o]), 
          (Subscript[kn, o])} = {(Subscript[Qx, o]), (Subscript[Qn, o]), (Subscript[kx, o]), 
            (Subscript[kn, o])} /. (Solve[{eq1, eq2, eq3, eq4}, {(Subscript[Qx, o]), (Subscript[Qn, o]), (Subscript[kx, o]), 
                  (Subscript[kn, o])}][[1]]); 
      eq5 := (((1 - (Subscript[θ, kxo])) (Subscript[Qx, o]))/(Subscript[lx, o])) == (((1 - (Subscript[θ, kno])) (Subscript[Pn, o]) (Subscript[Qn, o]))/(Subscript[ln, o])); (Subscript[lx, o]) =. ; 
      (Subscript[lx, o]) = (Subscript[lx, o]) /. (Solve[{eq5}, {(Subscript[lx, o])}][[1]]); (Subscript[k, 0]) = 1; 
      (Subscript[k, 1]) =. ; {(Subscript[k, 1])} = {(Subscript[k, 1])}  /.  
          FindRoot[{-(Subscript[V, 1])[(Subscript[Pn, o]), (Subscript[Ex, o])]/(Subscript[V, 2])[(Subscript[Pn, o]), (Subscript[Ex, o])] == 
                (Subscript[Cn, o])}, {(Subscript[k, 1]), (10^(-3))}, MaxIterations -> 1000, 
            AccuracyGoal -> $MinPrecision, 
            WorkingPrecision -> digits]; 
      (Subscript[a, 0]) =. ; 
{(Subscript[a, 0])} = {(Subscript[a, 0])}  /.  
          FindRoot[{ (Subscript[ϕ, 2])[(Subscript[Pn, o]), (Subscript[m, o])] ==  
                (Subscript[V, 2])[(Subscript[Pn, o]), (Subscript[Ex, 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[kx, o]), (Subscript[lx, o]), 
                  (Subscript[z, o])]  ==  (ρ + δ) 
                  (Subscript[Pk, o]), (Subscript[θ, kxo]) == (((ρ + δ) (Subscript[Pk, o]) (Subscript[kx, o]))/F)[(Subscript[kx, o]), (Subscript[lx, o]), (Subscript[z, o])]}, {(Subscript[a, 1]), 1}, {(Subscript[a, 2]), 1}, 
            MaxIterations -> 1000, AccuracyGoal -> digits, 
            WorkingPrecision -> digits]; {(Subscript[a, 3]), 
          (Subscript[a, 4])} =. ; 
{(Subscript[a, 3]), (Subscript[a, 4])} = {(Subscript[a, 3]), (Subscript[a, 4])} /. 
            FindRoot[{(Subscript[Pn, o]) (Subscript[G, 2])[(Subscript[kn, o]), (Subscript[ln, o]), (Subscript[z, o])]  == 
                  (Subscript[F, 2])[(Subscript[kx, o]), (Subscript[lx, o]), (Subscript[z, o])] , 
                1 - (Subscript[θ, kno]) == (((Subscript[G, 2])[(Subscript[kn, o]), (Subscript[ln, o]), (Subscript[z, o])] (Subscript[ln, o]))/G)[(Subscript[kn, o]), (Subscript[ln, o]), (Subscript[z, o])]}, {(Subscript[a, 3]), 1}, {(Subscript[a, 4]), 1}, 
              MaxIterations -> 1000, AccuracyGoal -> digits, 
              WorkingPrecision -> digits];

      {{(Subscript[m, o]), (Subscript[Ex, o]), (Subscript[Pk, o]), (Subscript[Pn, o]), (Subscript[kx, o]), (Subscript[ln, o]), (Subscript[lx, o]), (Subscript[Y, o]), (Subscript[Qn, o]), (Subscript[Qx, o]), 
              (Subscript[b, o]), g, (Subscript[z, o]), (Subscript[π, o]), (Subscript[Ex, o])/(Subscript[m, o]), (Subscript[m, o])/(Subscript[Y, o])}} // N // 
        MatrixForm

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

Solving for the Global Nonlinear Saddlepath

z0=(Subscript[z, o]);Pk[Pn_]=(Pn^α);Cn[Pn_, Ex_] = -(Subscript[V, 1])[Pn, Ex]/(Subscript[V, 2])[Pn, Ex];Ex[Pn_, ϕ1_]=Ex/.Solve[(Subscript[V, 2])[Pn, Ex]==ϕ1,Ex][[4]];
ix[Pn_, ϕ1_, ϕ2_, kx_] = 
              ix /.  Solve[ϕ2 == ϕ1 Pk[
                          Pn] (1 + ψ'[
                              ix/kx]), {ix}][[1]];
lx[Pn_, ϕ1_, ϕ2_, kx_, z_]=lx/. Solve[Cn[Pn, 
                            Ex[Pn, ϕ1]] + ((α Pk[Pn] (ix[Pn, ϕ1, ϕ2, kx] + ψ[ix[Pn, ϕ1, ϕ2, kx]/kx] kx))/Pn)==G[(Subscript[kn, o]),lbar - lx,z0],{lx}][[1]];
Pn[ϕ1_?NumberQ, ϕ2_?NumberQ, kx_?NumberQ, 
                z_?NumberQ]:=Block[{Pn},{Pn}={Pn}/.FindRoot[{Pn (Subscript[G, 2])[(Subscript[kn, o]),lbar - 
                                        lx[Pn, ϕ1, ϕ2, kx, z0],z0]==(Subscript[F, 2])[kx,lx[Pn,ϕ1,ϕ2,kx,(Subscript[z, o])],z]},{Pn, Pnguess},MaxIterations -> 1000,WorkingPrecision -> digits];Pnguess = Pn;Pn];
endovarsubs={Pk -> Pk[Pn],Ex -> Ex[Pn, ϕ1],ix -> ix[Pn, ϕ1, ϕ2, kx],lx->lx[Pn,ϕ1,ϕ2,kx,(Subscript[z, o])]};

Pndot:=(Dt[Pn]/.Solve[{Dt[Pn (Subscript[G, 2])[(Subscript[kn, o]),lbar - lx,z0]==(Subscript[F, 2])[kx, lx, z]],
Dt[Cn[Pn, 
                                      Ex] + ((α Pk[Pn] (ix + ψ[ix/kx] kx))/Pn)==G[(Subscript[kn, o]),lbar - lx,z0]],Dt[ϕ2 == ϕ1 Pk[
                                      Pn] (1 + ψ'[ix/kx])],Dt[(Subscript[V, 2])[Pn, Ex]==ϕ1]},{Dt[Pn], Dt[lx], Dt[ix], Dt[Ex]}][[1]])/.{Dt[ϕ1] -> ϕ1dot, 
                  Dt[ϕ2] -> ϕ2dot, Dt[kx] -> kxdot};
ϕ1dot := 
            0; ϕ2dot :=  (ρ + δ) ϕ2 - ϕ1 ((Subscript[F, 1])[kx, lx, z] + Pk (ψ'[ix/kx] ix/kx - ψ[ix/kx]));
bdot:= Pn G[(Subscript[kn, o]),lbar - lx,z0]+F[kx, lx, z]+ρ b-Ex-Pk (ix + ψ[ix/kx] kx);kxdot := ix - δ kx;z[t_] := (Subscript[z, temp])+((Subscript[z, o])-(Subscript[z, temp]))UnitStep[t - ts];
system:={Pndot,ϕ1dot,ϕ2dot,bdot,kxdot}/.endovarsubs;

(* 
          Set up system and initial conditions for the ODE solver *) 
sys:=Thread[{Pn'[t],ϕ1'[t],ϕ2'[t],b'[t],kx'[t]}== (system/.{Pn->Pn[t],ϕ1->ϕ1[t],ϕ2->ϕ2[t],b->b[t],kx->kx[t],z->z[t]}),List];
ini:=Thread[{Pn[0],ϕ1[0],ϕ2[0],b[0],kx[0]}=={Pnin,ϕ1in,ϕ2in,bin,kxin} ,List];
ini2:=Thread[{Pn[ts],ϕ1[ts],ϕ2[ts],b[ts],kx[ts]}=={Pnin2,ϕ1in2,ϕ2in2,bin2,kxin2} ,List];
  
box:=Ex[Pn[t], ϕ1[t]]>((Subscript[Ex, o])/5)||Ex[Pn[t], ϕ1[t]]<((2 (Subscript[Ex, o]))/1)|| ix[Pn[t], ϕ1[t], ϕ2[t], kx[t]]>((δ (Subscript[kx, o]))/5)||ix[Pn[t], ϕ1[t], ϕ2[t], kx[t]]<((2 δ (Subscript[kx, o]))/1)|| lx[Pn[t], ϕ1[t], ϕ2[t], kx[t], z[t]]>(1/5)||lx[Pn[t], ϕ1[t], ϕ2[t], kx[t], z[t]]<(8/10)||b[t]>-3 (Subscript[Y, o])||b[t]<3  (Subscript[Y, o])|| kx[t]>((Subscript[kx, o])/5)||kx[t]<(4/2) (Subscript[kx, o]); 
  
sysdiff[t1_, t2_]:={Max[Abs[Table[(Ex[Pn[t], ϕ1[t]]/.If[t < ts, gminsol1, gminsol2])-(Ex[Pn[t], ϕ1[t]]/.If[t < ts, gminsol11, gminsol21]),{t, t1, t2, 1/10}]]],Max[Abs[Table[(ix[Pn[t], ϕ1[t], ϕ2[t], 
                                  kx[t]]/.If[t < ts, gminsol1, gminsol2])-(ix[Pn[t], ϕ1[t], ϕ2[t], 
                                  kx[t]]/.If[t < ts, gminsol11, gminsol21]),{t, t1, t2, 1/10}]]],
Max[Abs[Table[(b[t]/.If[t < ts, gminsol1, gminsol2])-(b[t]/.If[t < ts, gminsol11, gminsol21]),{t, t1, t2, 1/10}]]],Max[Abs[Table[(kx[t]/.If[t < ts, gminsol1, gminsol2])-(kx[t]/.If[t < ts, gminsol11, gminsol21]),{t, t1, t2, 1/10}]]]};

distance[sol1_, sol2_, tsol_]:=
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) > 
                      Max[tol, gmindis] (10^(-5)),dislast = dis;traj=Table[Norm[Flatten[{Ex[Pn[t], ϕ1[t]] - 
                                        (Subscript[Ex, f]),ix[
                                        Pn[t], ϕ1[t], ϕ2[t], 
                                        kx[t]]-δ (Subscript[kx, f]),b[t]-(Subscript[b, f]),kx[t]-(Subscript[kx, f]),1000(z[t]-(Subscript[z, f]))}/.If[t < ts, sol1, sol2]],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}];

\(difftmin := 
            Print[\*""Difference in values of variables between two paths with (Subscript[t, min])s differing by 5 years = "", sysdiff[0, T] // N];\)
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] // N];
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."];
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 MAIN system!";

CIRCLE SEARCH

(* 
          Set up for solving for the global nonlinear saddle path *) 
Clear[compdis]; Clear[compdis2]; Off[Precision::precsm];
compdis[p_, epsilon_]:=Block[{ϕ1in,Pnin,ϕ2in,bin,kxin,zin,ϕ1min1,Pnmin1,ϕ2min1,ϕ1min2,Pnmin2,ϕ2min2,mindis,j1,j2,firstiter,sol,tsol,dis,tdis,tmin,jmin,minsol1,minsol2,tmin2,jmin2,minsol12,minsol22,mindis2},mindis = mindis2 = (10^7);ntraj = ntrajini;firstiter = 1;(Subscript[kx, f]) = (Subscript[kx, o]);(Subscript[lx, f]) = (Subscript[lx, o]);
 (* 
                      Pre - process the system for faster execution *) 
{Pnin,ϕ1in,ϕ2in,bin,kxin}:={(Subscript[Pn, o]),(Subscript[ϕ1, o]),(Subscript[ϕ2, o]),(Subscript[b, o]),(Subscript[kx, o])};
state=First[NDSolve`ProcessEquations[{sys, ini},{Pn, ϕ1, ϕ2, b, kx},t,WorkingPrecision -> digits ,MaxSteps -> 100,\(Method -> {"EventLocator", \ "Event" -> box, "EventAction" :> \ 
                                Throw[
                                  Null, \ "StopIntegration"], "EventLocationMethod" -> "StepBegin"}\)]];solve:=(theta = ((2  π j)/ntraj);{ϕ1in,ϕ2in}=p + epsilon {Sin[theta], Cos[theta]};{bin,kxin}={(Subscript[b, o]),(Subscript[kx, o])};
  
Pnguess = (Subscript[Pn, o]);Pnin=Pn[ϕ1in, ϕ2in, kxin, z[0]];
  
ϕ(Subscript[1, f]) = ϕ1in;{kxguess, lxguess} = {(Subscript[kx, f]), (Subscript[lx, f])};{(Subscript[kx, f]), (Subscript[lx, f])} =. ;(Subscript[z, f]) = (Subscript[z, temp]) + ((Subscript[z, o]) - (Subscript[z, temp]));(Subscript[Pn, f]) = 
                            (Subscript[F, 2])[(Subscript[kx, f]), (Subscript[lx, f]), (Subscript[z, f])]/(Subscript[G, 2])[(Subscript[kn, o]), lbar - (Subscript[lx, f]), (Subscript[z, f])];ϕ(Subscript[2, f]) = ϕ(Subscript[1, f]) Pk[(Subscript[Pn, f])];(Subscript[Ex, f]) = Ex[(Subscript[Pn, f]), ϕ(Subscript[1, f])];
{(Subscript[kx, f]), (Subscript[lx, f])} = {(Subscript[kx, f]), (Subscript[lx, f])} /. 
                              FindRoot[{Cn[(Subscript[Pn, f]), 
                                        Ex[(Subscript[Pn, f]), ϕ(Subscript[1, f])]] + α Pk[
                                        (Subscript[Pn, f])] δ (Subscript[kx, f])/
                                        (Subscript[Pn, f]) == 
                                    G[(Subscript[kn, o]), lbar - (Subscript[lx, f]), 
                                      (Subscript[z, f])], (ρ + δ) 
                                      Pk[(Subscript[Pn, f])] == 
                                    (Subscript[F, 1])[(Subscript[kx, f]), (Subscript[lx, f]), (Subscript[z, f])]}, {(Subscript[kx, f]), 
                                  kxguess}, {(Subscript[lx, f]), lxguess}, 
                                MaxIterations -> 1000, 
                                WorkingPrecision -> digits];
(Subscript[b, f]) = (((Subscript[Ex, f]) + Pk[(Subscript[Pn, f])] δ (Subscript[kx, f]) - F[(Subscript[kx, f]), (Subscript[lx, f]), (Subscript[z, f])] - (Subscript[Pn, f]) G[(Subscript[kn, o]), lbar - (Subscript[lx, f]), (Subscript[z, f])])/ρ);
{Pn[0],ϕ1[0],ϕ2[0],b[0],kx[0]}={Pnin,ϕ1in,ϕ2in,bin,kxin} ;If[(box /. t -> 0) ≠  
                              True,Print["Out of Box."];{dis, tdis} = {(10^7), 100};Goto[skip]];{Pn[0],ϕ1[0],ϕ2[0],b[0],kx[0]}=.;state1  = 
                          SetPrecision[
                            First[ NDSolve`Reinitialize[state,  ini]], 
                            precision];NDSolve`Iterate[
                          state1, { 0, ts - (10^(-12))}];sol1 = NDSolve`ProcessSolutions[state1];
If[
                          sol1[[1, 2, 1, 1, 2]] < ts - (10^(-12)), 
                          tsol = sol1[[1, 2, 1, 1, 2]]; {dis, tdis} = 
                            distance[sol1, sol2, tsol];];
If[sol1[[1, 2, 1, 1, 2]] ==  
                              ts - (10^(-12)),{ϕ1in2,ϕ2in2,bin2,kxin2}={ϕ1[ts-(10^(-12))]/.sol1,ϕ2[ts-(10^(-12))]/.sol1,b[ts-(10^(-12))]/.sol1,kx[ts-(10^(-12))]/.sol1};
  
Pnguess = 
                                  Pn[ts - (10^(-12))] /. sol1;Pnin2=Pn[ϕ1[ts-(10^(-12))]/.sol1,ϕ2[ts - (10^(-12))] /. 
                                        sol1,kx[ts - (10^(-12))] /. sol1,z[ts]];state2  = 
                                SetPrecision[ 
                                  First[NDSolve`Reinitialize[state,  ini2]], 
                                  precision];NDSolve`Iterate[state2, { ts, Tmax}];sol2 = 
                                NDSolve`ProcessSolutions[state2];
{dis, tdis} = 
                                distance[sol1, sol2, 
                                  sol2[[1, 2, 1, 1, 2]]];];Label[skip];);
If[epsilon == 0, solve; 
                    mindis2 = dis; tmin2 = tdis; Goto[end1]];
For[j = 1,j ≤  ntraj  ,j++,If[trajtrack == 
                              0 || (trajtrack == 
                                  1  && (j == minpos[[1]] || 
                                    j == minpos[[2]])),solve;jdis[j] = dis;jtdis[j] = tdis;dissol1[j] = sol1;dissol2[j] = sol2;jp[j]={ϕ1in,ϕ2in},jdis[j] = (10^7)]];tab = 
                    Table[jdis[Mod[ii - 1, ntraj] + 1] - 
                        jdis[Mod[ii, ntraj] + 1], {ii, 1, ntraj}];minpos = 
                    Flatten[
                      Position[
                        Table[tab[[Mod[ii - 2, ntraj] + 1]] > 0  && 
                            tab[[Mod[ii - 1, ntraj] + 1]] < 0, {ii, 1, 
                            ntraj}], True]];
For[ii = 1, 
                    ii ≤ Length[minpos], ii++, 
                    Print["ii = ", 
                      ii, "   {mindis, tdis, jmin, ntraj} = ", {jdis[
                          minpos[[ii]]], jtdis[minpos[[ii]]] // N, 
                        minpos[[ii]], ntraj}]];If[
                    Length[minpos] == 2  &&  iflag == 1, 
                    trajtrack = 1];firstiter = 0;

For[ii = 1,ii ≤ Length[minpos],ii++,{mindis, tmin, jmin, ntraj, minsol1, 
                            minsol2} = {jdis[minpos[[ii]]], 
                            jtdis[minpos[[ii]]], minpos[[ii]], 
                            ntrajini, dissol1[minpos[[ii]]], 
                            dissol2[minpos[[ii]]]};{ϕ1min1,ϕ2min1}=jp[minpos[[ii]]];{j1, j2} = 2 {jmin - 1, jmin + 1};{ntraj, jmin} = 2 {ntraj, jmin};tminmin = 0;Label[begin];
For[j = j1 + 1,j ≤  j2  ,j += 2,solve;If[ Abs[mindis - dis] ≤ 
                                    Max[(tol/100), (gmindis/100)],Print[mindis]; Goto[end]; ];If[dis ≤ mindis,mindis = dis;tmin = tdis;jmin = j;minsol1 = sol1;minsol2 = sol2;{ϕ1min1,ϕ2min1}={ϕ1in,ϕ2in};Break[];];];
{j1, j2} = 
                          2 {jmin - 1, jmin + 1};{ntraj, jmin} = 2 {ntraj, jmin}; Goto[begin];Label[end];
If[mindis < mindis2,mindis2 = mindis;tmin2 = tmin;jmin2 = jmin;minsol12 = minsol1;minsol22 = minsol2;{ϕ1min2,ϕ2min2}={ϕ1min1,ϕ2min1}];
If[
                          mindis2 < gmindis  &&  
                            gmindis ≠  (10^7), Break[]];];
Print["This Iteration:    \
Minimum Distance = ", 
                    mindis2 // 
                      N, "    No. of Subdivisions of ϵ-ring > ", 
                    jmin2];Label[end1];{{mindis2, tmin2, minsol12, minsol22},{ϕ1min2,ϕ2min2}}];

compdis2[p_, ϵ_]:=Block[{},iflag = 0;If[
                    firstsim == 
                      1, {pmin, epsilon} = {p, ϵ}, {pmin, 
                        epsilon} = (<< solution2)[[{5, 6}]]; 
                    iflag = 1];trajtrack = 0;rit = 1;rinfo = {};Label[begin];Print["\nBEGINNING ITERATION ",rit," WITH CENTER OF ϵ-ring at ",SetPrecision[pmin, precision]];{{rmin, rtmin, rminsol1, rminsol2},{ϕ1rmin,ϕ2rmin}}=compdis[pmin, epsilon];If[rmin ≤ gmindis,gmindis = rmin;gtmin = rtmin;gminsol1 = rminsol1;gminsol2 = rminsol2;pmin={ϕ1rmin,ϕ2rmin};If[firstsim == 1,epsilon=epsilon/1],epsilon = epsilon/2];Print["AT THE END OF ITERATION = ",rit,":    CURRENT MINIMUM DISTANCE = ",gmindis // N,"     TIME OF CLOSEST APPROACH = ",gtmin];
If[firstsim == 1,If[gtmin > T + 5  &&  iflag == 0, 
                          iflag = 
                            1;   ; {gminsol1, gminsol2, 
                              gtmin, gmindis, pmin, 5  epsilon} >> solution1; 
                          Print[""]; Print[step2];];
If[gtmin > 
                                Max[T + 5 , (<< solution1)[[3]] + 
                                    5] &&  iflag == 1 ,{gminsol11, gminsol21, 
                                  gtmin1} = (<< solution1)[[{1, 2, 
                                      3}]];If[Max[sysdiff[0, T]] > tol,Print[""]; difftmin; 
                                  Print[
                                    notsatistmin]; {gminsol1, gminsol2, 
                                      gtmin, gmindis, pmin, 5  epsilon} >> 
                                    solution1,Print[""];difftmin;Print[satistmin];{gminsol1, gminsol2, gtmin, 
                                        gmindis, pmin, 5  epsilon} >> 
                                      solution1;Print[""];Print[step3];Goto[end]]],
If[gtmin > (<< solution2)[[3]] &&  
                              gmindis < (<< solution2)[[4]],{gminsol1, gminsol2, gtmin, gmindis, 
                                  pmin, 5  epsilon} >> solution1;Goto[end]];];
AppendTo[
                    rinfo, {gtmin, gmindis, pmin, 5  epsilon}];If[Length[rinfo]>15 && (((Part[rinfo, Length[rinfo] - 10][[2]] - Last[rinfo][[2]])/(Part[rinfo, Length[rinfo] - 10]))[[2]]) < (10^(-1))&&(((Part[rinfo, Length[rinfo] - 10][[1]] - Last[rinfo][[1]])/(Part[rinfo, Length[rinfo] - 10]))[[1]]) < (10^(-1)),Print["You may be losing precision in calculations."];gtmin = (10^100);epsilon = Part[rinfo, 
                              Length[rinfo] - 10][[4]];pmin = Part[rinfo, 
                              Length[rinfo] - 10][[3]];Goto[end];];rit++;Goto[begin];Label[end];{gminsol1, gminsol2, gtmin, 
                    gmindis, pmin, 5  epsilon}];

(* INITIALIZATION CELL - Initialization for finding the required trajectory. *)
{(Subscript[ϕ1, o]),(Subscript[ϕ2, o])}={(Subscript[V, 2])[(Subscript[Pn, o]), (Subscript[Ex, o])], (Subscript[V, 2])[(Subscript[Pn, o]), (Subscript[Ex, o])] Pk[(Subscript[Pn, o])]};          
T = 20;                   
ts=5;                   ;
tol=(10^(-4));           ;
epsilon= (((10^(-0)))/4)(((((Subscript[ϕ1, o])^2)+((Subscript[ϕ2, o])^2)))^((1/2)));                                                  ;
pini={(9/10)(Subscript[ϕ1, o]),(11/10)(Subscript[ϕ2, o])};

Timing[precision=16;firstsim = 1;Tmax = 200;ϵmax = epsilon;pinit = pini;gmindis = 1;
While[True,localmin = 0;Label[begin];$MinPrecision = digits = precision;ntrajini = 20;firstiter = 1;If[
                firstsim == 
                  0, {pini, epsilon} = (<< solution2)[[{5, 6}]]];If[localmin == 1, 
                pini = pinit = pinit + ϵmax/10 {1/2, 1/2}];
  
{gmindis, gtmin} = 
                Part[Part[compdis[pini, 0], 1], {1, 2}];Print["AT THE CENTER OF THE INITIAL RING WITH CURRENT PRECISION: ","    MINIMUM DISTANCE = ",gmindis,"    TIME OF CLOSEST APPROACH = ",gtmin];
  ;
{gminsol1, gminsol2, gtmin, gmindis, pmin, eps} = 
                compdis2[pini, epsilon];If[gtmin == (10^100),precision = precision + 4;epsilon = SetPrecision[eps, precision];pini = SetPrecision[pmin, precision];If[precision < 80,increaseprec,Print["Either the system is highly ill-conditioned or you might have encountered a local minima."];Print["Perturing the center of the initial ring and resetting 'precision' to its initial value."];precision = 16;localmin = 1];Goto[begin]];
   
If[firstsim  ==  1,firstsim = 0; 
                  precision = precision + 4; 
                  increaseprec; {gminsol1, gminsol2, gtmin, gmindis, pmin, 
                      eps} >> solution2,
{gminsol11, gminsol21, 
                        gtmin1} = (<< solution2)[[{1, 2, 3}]];If[Max[sysdiff[0, T]] > tol,Print[""]; diffprec; Print[notsatisprec]; 
                        precision = precision + 4; 
                        increaseprec; {gminsol1, gminsol2, gtmin, gmindis, 
                            pmin, eps} >> solution2,Print[""];diffprec;Print[satisprec];Print[""];Print[congratulations];Goto[end]]];];   ;
Label[end];For[j = T + 1,j ≤ IntegerPart[Min[gtmin1, gtmin]],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[gtmin1, gtmin]] + 1,Print["The solution is accurate upto ",j - 1," years. Continue with next cell."]]]

LINE SEARCH for m

{gminsol1, gminsol2, 
              gtmin} = (<< solution2)[[{1, 2, 3}]];
ϕ(Subscript[1, f]) = ϕ1[0] /. gminsol1; 
          (Subscript[Pn, f]) = (Subscript[F, 2])[(Subscript[kx, f]), (Subscript[lx, f]), (Subscript[z, f])]/(Subscript[G, 2])[(Subscript[kn, o]), lbar - (Subscript[lx, f]), (Subscript[z, f])]; 
          (Subscript[ln, f]) = lbar - (Subscript[lx, f]); (Subscript[Qn, f]) = G[(Subscript[kn, o]), lbar - (Subscript[lx, f]), (Subscript[z, f])]; 
          (Subscript[Qx, f]) = F[(Subscript[kx, f]), (Subscript[lx, f]), (Subscript[z, f])]; (Subscript[Y, f]) = (Subscript[Pn, f]) (Subscript[Qn, f]) + (Subscript[Qx, f]); 
          (Subscript[Pk, f]) = Pk[(Subscript[Pn, f])]; {(Subscript[kx, f]), (Subscript[lx, f])} =. ;
{(Subscript[kx, f]), (Subscript[lx, f])} = {(Subscript[kx, f]), (Subscript[lx, f])} /. 
                FindRoot[{Cn[(Subscript[Pn, f]), 
                          Ex[(Subscript[Pn, f]), ϕ(Subscript[1, f])]] + α Pk[
                            (Subscript[Pn, f])] δ (Subscript[kx, f])/(Subscript[Pn, f]) == 
                      G[(Subscript[kn, o]), lbar - (Subscript[lx, f]), (Subscript[z, f])], (ρ + δ) 
                        Pk[(Subscript[Pn, f])] == (Subscript[F, 1])[(Subscript[kx, f]), (Subscript[lx, f]), (Subscript[z, f])]}, {(Subscript[kx, f]), 
                    (Subscript[kx, o])}, {(Subscript[lx, f]), (Subscript[lx, o])}, MaxIterations -> 1000, 
                  WorkingPrecision -> digits];
(Subscript[Ex, f]) = 
            Ex[(Subscript[Pn, f]), ϕ(Subscript[1, f])]; ϕ(Subscript[2, f]) = ϕ(Subscript[1, f]) Pk[(Subscript[Pn, f])]; 
          (Subscript[b, f]) = (((Subscript[Ex, f]) + Pk[(Subscript[Pn, f])] δ (Subscript[kx, f]) - F[(Subscript[kx, f]), (Subscript[lx, f]), (Subscript[z, f])] - (Subscript[Pn, f]) G[(Subscript[kn, o]), lbar - (Subscript[lx, f]), (Subscript[z, f])])/ρ);
{(Subscript[m, f]), (Subscript[π, f])} =. ;
{(Subscript[m, f]), (Subscript[π, f])} = {(Subscript[m, f]), (Subscript[π, f])}  /.  
                FindRoot[{ (Subscript[ϕ, 2])[(Subscript[Pn, f]), (Subscript[m, f])] ==  
                      (Subscript[V, 2])[(Subscript[Pn, f]), (Subscript[Ex, 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];
\(\({{"m", "Ex", "Pk", "Pn", "kn", "kx", "ln", "lx", "Y", "Qn", "Qx", "b", "g", "z", "π", "Ex/m", "m/Y"}, {m\_o, Ex\_o, Pk\_o, Pn\_o, 
                  kn\_o, kx\_o, ln\_o, lx\_o, Y\_o, Qn\_o, Qx\_o, b\_o, g, 
                  z\_o, π\_o, Ex\_o/m\_o, m\_o/Y\_o}, {m\_f, Ex\_f, 
                  Pk\_f, Pn\_f, kn\_o, kx\_f, ln\_f, lx\_f, Y\_f, Qn\_f, 
                  Qx\_f, b\_f, g, z\_f, π\_f, Ex\_f/m\_f, m\_f/Y\_f}} // 
              N\) // MatrixForm\)

Clear[χ]; χ[m_?NumberQ, 
            t_?NumberQ] := (( (Subscript[ϕ, 2])[Pn[t], m])/( ϕ1[t])) - ρ /. If[t < ts, gminsol1, gminsol2]

distancem[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, lmindis] (10^(-5)),dislast = dis;traj=Table[Norm[Flatten[{m[t]-(Subscript[m, 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 = 10];{dis, tdis}];
  
sysdiffm[t1_, t2_] := 
            Max[Abs[Table[(m[t] /. lminsol2) - (m[t] /. 
                        lminsol21), {t, t1, t2, 1/10}]]];

\(difftminm := 
            Print[\*""Difference in values of variables between two paths with (Subscript[t, min])s differing by 5 years = "", sysdiffm[ts, T] // N];\)
notsatistminm="The (Subscript[(t), min]) is NOT large enough. The variables changed by greater than\! 'tol' due to latest increase in (Subscript[(t), min]).";
satistminm= "THE ACCURACY CHECK PASSED SUCCESSFULLY. The variables changed by less\! than\! 'tol' due to latest increase in (Subscript[(t), min]).";
congratulationsm= "CONGRATULATIONS! You are COMPLETELY done! Hope to see you in another model.";

(* 
          Set up for solving for the global nonlinear saddle path *) 
Clear[compdism];
compdism[p_, ϵ_]:=Block[{pmin, epsilon},  
iflag = 0;lit = 1;{pmin, epsilon} = {p, ϵ};Label[begin];
If[epsilon ≠ 0,Print["\nBEGINNING ITERATION ",lit," WITH LENGTH and CENTER OF ϵ-segment = ",SetPrecision[{epsilon, pmin}, 
                              precision]];];
For[j = -1,j ≤  1  ,j = j + 2,min2=pmin + j epsilon;sol = First[NDSolve[SetPrecision[{m'[t] == 
                                        g - χ[m[t], t] m[t] ,m[ts]==min2},precision],{m},{t, ts, Tmax},Method->{"EventLocator", "Event"->m[t] > Min[(Subscript[m, o]), (Subscript[m, f])]/2 || 
                                        m[t] < 2  Max[(Subscript[m, o]), (Subscript[m, f])],\("EventAction" :> \ Throw[Null, \ "StopIntegration"]\),\("EventLocationMethod" -> "StepBegin"\)},WorkingPrecision -> digits ,MaxSteps -> 300]];{dis, tdis} = 
                          distancem[sol, sol[[1, 2, 1, 1, 2]]];jdis[j] = dis;jtdis[j] = tdis;dissol[j] = sol;jp[j]=min2;jm[j]=m[sol[[1, 2, 1, 1, 2]]]/.sol;];
If[jm[1]>(Subscript[m, f]) && jm[-1]>(Subscript[m, f]),{lmindis, ltmin, lminsol2, 
                          pmin} = {jdis[-1], jtdis[-1], 
                          dissol[-1], jp[-1]},If[jm[1]>(Subscript[m, f]),epsilon = epsilon/2,{lmindis, ltmin, lminsol2, pmin} = {jdis[1], 
                              jtdis[1], dissol[1], jp[1]}]];

If[epsilon == 0,lmindis = jdis[-1];ltmin = jtdis[-1];Goto[end1]];

Print["AT THE END OF ITERATION = ",lit,"     TIME OF CLOSEST APPROACH = ",ltmin,"     MINIMUM DISTANCE = ",lmindis];
If[
                    ltmin > T + 5  &&  iflag == 0, 
                    iflag = 1; {lminsol2, ltmin, lmindis, pmin} >> 
                      solutionm1; Print[""]; Print[step2];];
If[ltmin > 
                          Max[T + 5 , (<< solutionm1)[[2]] + 5] &&  
                        iflag == 1 ,{lminsol21, 
                            ltmin1} = (<< solutionm1)[[{1, 2}]];If[sysdiffm[ts, T] > tol,Print[""]; difftminm; 
                            Print[notsatistminm]; {lminsol2, ltmin, lmindis, 
                                pmin} >> solutionm1,Print[""];difftminm;Print[satistminm];{lminsol2, ltmin, lmindis, pmin} >> 
                                solutionm1;Goto[end]]];lit++;Goto[begin];Label[end];
  ;
 (* 
                      Solve for the value of m at ts - . In many cases it may not jump at ts as in this example *) 
min=m[ts]/.lminsol2;
lminsol1 = First[NDSolve[SetPrecision[{m'[t] == 
                                      g - χ[m[t], t] m[t] ,m[ts - (10^(-12))]==min},precision],{m},{t, ts - (10^(-12)), 0},WorkingPrecision -> digits ,MaxSteps -> 300]];{lminsol1, lminsol2, ltmin, lmindis, pmin} >> 
                    solutionm1;
  ;
For[j = T + 1,j ≤ IntegerPart[Min[ltmin1, ltmin]],j++,
If[Max[sysdiffm[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[ltmin1, ltmin]] + 1,Print["The solution is accurate upto ",j - 1," years. Continue with next cell."]];Label[end1];
{lminsol1, lminsol2, ltmin, 
                    lmindis, pmin}];

Timing[Tmax = gtmin;$MinPrecision = digits = precision;lmindis = 1;
  ;{lmindis, ltmin}=Part[compdism[(Subscript[m, o]),0],{4, 3}];Print[{lmindis, ltmin}];
Print["AT THE CENTER OF THE INITIAL SEGMENT: ","    MINIMUM DISTANCE = ",lmindis,"    TIME OF CLOSEST APPROACH = ",ltmin];
  ;{lminsol1, lminsol2, ltmin, lmindis, pmin}=compdism[(Subscript[m, o]),((Subscript[m, o])/4)];Print[congratulationsm]; ]

RESULTS

{temp1, temp2} = << fig13; T = 20; (Subscript[z, f]) = 1;

{gminsol1, gminsol2, 
                gtmin} =   
                temp2[[{1, 2, 3}]];;
ϕ(Subscript[1, f]) = ϕ1[0] /. gminsol1; 
          (Subscript[Pn, f]) = (Subscript[F, 2])[(Subscript[kx, f]), (Subscript[lx, f]), (Subscript[z, f])]/(Subscript[G, 2])[(Subscript[kn, o]), lbar - (Subscript[lx, f]), (Subscript[z, f])]; 
          (Subscript[ln, f]) = lbar - (Subscript[lx, f]); (Subscript[Qn, f]) = G[(Subscript[kn, o]), lbar - (Subscript[lx, f]), (Subscript[z, f])]; 
          (Subscript[Qx, f]) = F[(Subscript[kx, f]), (Subscript[lx, f]), (Subscript[z, f])]; (Subscript[Y, f]) = (Subscript[Pn, f]) (Subscript[Qn, f]) + (Subscript[Qx, f]); 
          (Subscript[Pk, f]) = Pk[(Subscript[Pn, f])]; {(Subscript[kx, f]), (Subscript[lx, f])} =. ;
{(Subscript[kx, f]), (Subscript[lx, f])} = {(Subscript[kx, f]), (Subscript[lx, f])} /. 
                FindRoot[{Cn[(Subscript[Pn, f]), 
                          Ex[(Subscript[Pn, f]), ϕ(Subscript[1, f])]] + α Pk[
                            (Subscript[Pn, f])] δ (Subscript[kx, f])/(Subscript[Pn, f]) == 
                      G[(Subscript[kn, o]), lbar - (Subscript[lx, f]), (Subscript[z, f])], (ρ + δ) 
                        Pk[(Subscript[Pn, f])] == (Subscript[F, 1])[(Subscript[kx, f]), (Subscript[lx, f]), (Subscript[z, f])]}, {(Subscript[kx, f]), 
                    (Subscript[kx, o])}, {(Subscript[lx, f]), (Subscript[lx, o])}, MaxIterations -> 1000, 
                  WorkingPrecision -> digits];
(Subscript[Ex, f]) = 
            Ex[(Subscript[Pn, f]), ϕ(Subscript[1, f])]; ϕ(Subscript[2, f]) = ϕ(Subscript[1, f]) Pk[(Subscript[Pn, f])]; 
          (Subscript[b, f]) = (((Subscript[Ex, f]) + Pk[(Subscript[Pn, f])] δ (Subscript[kx, f]) - F[(Subscript[kx, f]), (Subscript[lx, f]), (Subscript[z, f])] - (Subscript[Pn, f]) G[(Subscript[kn, o]), lbar - (Subscript[lx, f]), (Subscript[z, f])])/ρ);
{(Subscript[m, f]), (Subscript[π, f])} =. ;
{(Subscript[m, f]), (Subscript[π, f])} = {(Subscript[m, f]), (Subscript[π, f])}  /.  
                FindRoot[{ (Subscript[ϕ, 2])[(Subscript[Pn, f]), (Subscript[m, f])] ==  
                      (Subscript[V, 2])[(Subscript[Pn, f]), (Subscript[Ex, 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];
\(\({{"m", "Ex", "Pk", "Pn", "kn", "kx", "ln", "lx", "Y", "Qn", "Qx", "b", "g", "z", "π", "Ex/m", "m/Y"}, {m\_o, Ex\_o, Pk\_o, Pn\_o, 
                  kn\_o, kx\_o, ln\_o, lx\_o, Y\_o, Qn\_o, Qx\_o, b\_o, g, 
                  z\_o, π\_o, Ex\_o/m\_o, m\_o/Y\_o}, {m\_f, Ex\_f, 
                  Pk\_f, Pn\_f, kn\_o, kx\_f, ln\_f, lx\_f, Y\_f, Qn\_f, 
                  Qx\_f, b\_f, g, z\_f, π\_f, Ex\_f/m\_f, m\_f/Y\_f}} // 
              N\) // MatrixForm\)

Clear[χ]; χ[m_?NumberQ, 
            t_?NumberQ] := (( (Subscript[ϕ, 2])[Pn[t], m])/( ϕ1[t])) - ρ /. If[t < ts, gminsol1, gminsol2]

{{lminsol1, lminsol2}, {gminsol1, 
            gminsol2}} =   {temp1[[{1, 2}]], 
          temp2[[{1, 2}]]};
tplot =T;
pm=Plot[{(Subscript[m, o]),m[t]/.If[t < ts, lminsol1, lminsol2],(Subscript[m, f])},{t, 0, tplot},PlotRange -> All,AxesLabel->{"",m},DisplayFunction -> Identity];
pPn=Plot[{(Subscript[Pn, o]),Pn[t]/.If[t < ts, gminsol1, gminsol2],(Subscript[Pn, f])},{t, 0, tplot},PlotRange -> All,AxesLabel->{"",Pn},DisplayFunction -> Identity];
pEx=Plot[{(Subscript[Ex, o]),Ex[Pn[t], ϕ1[t]]/.If[t < ts, gminsol1, gminsol2],(Subscript[Ex, f])},{t, 0, tplot},PlotRange -> All,AxesLabel->{"",Ex},DisplayFunction -> Identity];
pix=Plot[{δ (Subscript[kx, o]),ix[Pn[t], ϕ1[t], ϕ2[t], kx[t]]/.If[t < ts, gminsol1, gminsol2],δ (Subscript[kx, f])},{t, 0, tplot},PlotRange -> All,AxesLabel->{"",ix},DisplayFunction -> Identity];
pb=Plot[{(Subscript[b, o]),b[t]/.If[t < ts, gminsol1, gminsol2],(Subscript[b, f])},{t, 0, tplot},PlotRange -> All,AxesLabel->{"",b},DisplayFunction -> Identity];
pkx=Plot[{(Subscript[kx, o]),kx[t]/.If[t < ts, gminsol1, gminsol2],(Subscript[kx, f])},{t, 0, tplot},PlotRange -> All,AxesLabel->{"",kx},DisplayFunction -> Identity];
pz=Plot[{(Subscript[z, o]),z[t]/.If[t < ts, gminsol1, gminsol2],(Subscript[z, f])},{t, 0, tplot},PlotRange -> All,AxesLabel->{"",z},DisplayFunction -> Identity];
pQx=Plot[{(Subscript[Qx, o]),F[kx[t], 
                      lx[Pn[t], ϕ1[t], ϕ2[t], kx[t], z0], z[t]]/.If[t < ts, gminsol1, gminsol2],(Subscript[Qx, f])},{t, 0, tplot},PlotRange -> All,AxesLabel->{"",Qx},DisplayFunction -> Identity];
pQn=Plot[{(Subscript[Qn, o]),G[(Subscript[kn, o]),lbar - 
                          lx[Pn[t], ϕ1[t], ϕ2[t], kx[t], z0],z0]/.If[t < ts, gminsol1, gminsol2],(Subscript[Qn, f])},{t, 0, tplot},PlotRange -> All,AxesLabel->{"",Qn},DisplayFunction -> Identity];
pY=Plot[{(Subscript[Y, o]),(Subscript[Pn, o])G[(Subscript[kn, o]),lbar - 
                              lx[Pn[t], ϕ1[t], ϕ2[t], kx[t], z0],z0]+F[kx[t], 
                        lx[Pn[t], ϕ1[t], ϕ2[t], kx[t], z0], z[t]]/.If[t < ts, gminsol1, gminsol2],(Subscript[Pn, o])(Subscript[Qn, f])+(Subscript[Qx, f])},{t, 0, tplot},PlotRange -> All,AxesLabel->{"",Y},DisplayFunction -> Identity];
plx=Plot[{(Subscript[lx, o]),lx[Pn[t], ϕ1[t], ϕ2[t], kx[t], z0]/.If[t < ts, gminsol1, gminsol2],(Subscript[lx, f])},{t, 0, tplot},PlotRange -> All,AxesLabel->{"",lx},DisplayFunction -> Identity];
pln=Plot[{lbar- (Subscript[lx, o]),lbar - 
                      lx[Pn[t], ϕ1[t], ϕ2[t], kx[t], z0]/.If[t < ts, gminsol1, gminsol2],lbar-(Subscript[lx, f])},{t, 0, tplot},PlotRange -> All,AxesLabel->{"",ln},DisplayFunction -> Identity];
pχ=Plot[{(Subscript[π, o]),χ[m[t], t]/.If[t < ts, lminsol1, lminsol2],(Subscript[π, f])},{t, 0, tplot},PlotRange -> All,AxesOrigin->{0,  .99  (Subscript[π, f])},AxesLabel->{"",χ},DisplayFunction -> Identity];pπ=Plot[{(Subscript[π, o]),χ[m[t], t]+((((((Subscript[k, 1])))^β) 
                                Pn[t]^(1 - β))/(((((Subscript[k, 0])))^β) + ((((Subscript[k, 1])))^β) Pn[t]^(1 - β)))((Pndot /. endovarsubs/.{Pn->Pn[t],ϕ1->ϕ1[t],ϕ2->ϕ2[t],b->b[t],kx->kx[t],z->z[t]})/(Pn[t]))/.If[t < ts, gminsol1, gminsol2]/.If[t < ts, lminsol1, lminsol2],(Subscript[π, f])},{t, 0, tplot},PlotRange -> All,AxesOrigin->{0,  .99  (Subscript[π, f])},AxesLabel->{"",π},DisplayFunction -> Identity];
Show[GraphicsArray[{{pm,pPn,pEx,pix},{pb,pkx,pz},{pQx,pQn,pY},{plx,pln,pπ,pχ}}]]

tplot =T;
pz=Plot[{If[t < 0,(Subscript[z, o]),z[t]/.If[t < ts, gminsol1, gminsol2]]},{t, -2, tplot},PlotRange -> { .95, 1.41},AxesOrigin -> {-2,  .95},AxesLabel -> {"t", "z"},PlotStyle -> {{Thickness[ .007], 
                  Dashing[{ .05, 0}]}},DisplayFunction -> Identity];pEx=Plot[{If[t < 0,((Subscript[Ex, o])/(c[(Subscript[Pn, o])])),((Ex[Pn[t], ϕ1[t]])/(c[Pn[t]]))/.If[t < ts, gminsol1, gminsol2]]},{t, -2, tplot},PlotRange -> { .245,  .275},AxesOrigin -> {-2,  .245},AxesLabel->{"t",(Subscript["E", 1])},PlotStyle -> {{Thickness[ .007], 
                  Dashing[{ .05, 0}]}},DisplayFunction -> Identity];pPn=Plot[{If[t < 0,(Subscript[Pn, o]),Pn[t]/.If[t < ts, gminsol1, gminsol2]]},{t, -2, tplot},PlotRange -> { .95, 1.41},AxesOrigin -> {-2,  .95},AxesLabel->{"t",(Subscript[P, n])},PlotStyle -> {{Thickness[ .007], 
                  Dashing[{ .05, 0}]}},DisplayFunction -> Identity];
pm=Plot[{If[t < 0,(Subscript[m, o]),m[t]/.If[t < ts, lminsol1, lminsol2]]},{t, -2, tplot},PlotRange -> { .08,  .115},AxesOrigin -> {-2,  .08},AxesLabel -> {"t", m},PlotStyle -> {{Thickness[ .007], 
                  Dashing[{ .05, 0}]}},DisplayFunction -> Identity];
pπ=Plot[{If[t < 0,(Subscript[π, o]),χ[m[t], 
                            t] + ((((((Subscript[k, 1])))^β) 
                                  Pn[t]^(1 - β))/(((((Subscript[k, 0])))^β) + ((((Subscript[k, 1])))^β) Pn[t]^(1 - β))) (1/Pn)[
                                  t] (Pndot /. 
                                  endovarsubs /. {Pn -> 
                                    Pn[t], ϕ1 -> ϕ1[
                                      t], ϕ2 -> ϕ2[t], 
                                  b -> b[t], kx -> kx[t], 
                                  z -> z[t]}) /. 
                        If[t < ts, gminsol1, gminsol2] /. 
                      If[t < ts, lminsol1, lminsol2]]},{t, -2, tplot},PlotRange -> { .36,  .53},AxesOrigin->{-2,.36},AxesLabel -> {"t", π},PlotStyle -> {{Thickness[ .007], 
                  Dashing[{ .05, 0}]}},DisplayFunction -> Identity]; pχ=Plot[{If[t < 0,(Subscript[π, o]),χ[m[t], t] /. 
                      If[t < ts, lminsol1, lminsol2]]},{t, -2, tplot},PlotRange -> { .36,  .53},AxesOrigin -> {-2,  .36},AxesLabel -> {"t", χ},PlotStyle -> {{Thickness[ .007], 
                  Dashing[{ .05, 0}]}},DisplayFunction -> Identity];
 plx=Plot[{If[t < 0,(Subscript[lx, o]),lx[Pn[t], ϕ1[t], ϕ2[t], kx[t], 
                          z[t]]/.If[t < ts, gminsol1, gminsol2]]},{t, -2, tplot},PlotRange -> { .45,  .55},AxesOrigin -> {-2,   .45},AxesLabel->{"t",(Subscript[L, x])},PlotStyle -> {{Thickness[ .007], 
                    Dashing[{ .05, 0}]}},DisplayFunction -> Identity];
pln=Plot[{If[t < 0,(Subscript[ln, o]),lbar - 
                        lx[Pn[t], ϕ1[t], ϕ2[t], kx[t], z[t]]/.If[t < ts, gminsol1, gminsol2]]},{t, -2, tplot},PlotRange -> { .45,  .55},AxesOrigin -> {-2,  .45},AxesLabel->{ "t",(Subscript[L, n])},PlotStyle -> {{Thickness[ .007], 
                  Dashing[{ .05, 0}]}},DisplayFunction -> Identity];  pQn=Plot[{If[t < 0,(Subscript[Qn, o]),G[(Subscript[kn, o]),lbar - 
                            lx[Pn[t], ϕ1[t], ϕ2[t], kx[t], z[t]],z[t]]/.If[t < ts, gminsol1, gminsol2]]},{t, -2, tplot},PlotRange -> { .50,  .81},AxesOrigin -> {-2,  .50},AxesLabel->{"t",(Subscript[Q, n])},PlotStyle -> {{Thickness[ .007], 
                  Dashing[{ .05, 0}]}},DisplayFunction -> Identity]; pix=Plot[{If[t < 0,δ (Subscript[kx, o]),ix[Pn[t], ϕ1[t], ϕ2[t], kx[t]]/.If[t < ts, gminsol1, gminsol2]]},{t, -2, tplot},PlotRange -> { .05,  .11},AxesOrigin -> {-2,  .05},AxesLabel->{"t",(Subscript["I", x])},PlotStyle -> {{Thickness[ .007], 
                  Dashing[{ .05, 0}]}},DisplayFunction -> Identity];pkx=Plot[{If[t < 0,(Subscript[kx, o]),kx[t]/.If[t < ts, gminsol1, gminsol2]]},{t, -2, tplot},PlotRange -> {1.65, 1.96},AxesOrigin -> {-2, 1.65},AxesLabel->{"t",(Subscript[K, x])},PlotStyle -> {{Thickness[ .007], 
                  Dashing[{ .05, 0}]}},DisplayFunction -> Identity];
pQx=Plot[{If[t < 0,(Subscript[Qx, o]),F[kx[t], 
                        lx[Pn[t], ϕ1[t], ϕ2[t], kx[t], z[t]], 
                        z[t]]/.If[t < ts, gminsol1, gminsol2]]},{t, -2, tplot},PlotRange -> { .45,  .85},AxesOrigin -> {-2,  .45},AxesLabel->{"t",(Subscript[Q, x])},PlotStyle -> {{Thickness[ .007], 
                  Dashing[{ .05, 0}]}},DisplayFunction -> Identity];pY=Plot[{If[t < 0,(Subscript[Pn, o]) (Subscript[Qn, o])+(Subscript[Qx, o]),(Subscript[Pn, o])G[(Subscript[kn, o]),lbar - 
                                lx[Pn[t], ϕ1[t], ϕ2[t], kx[t], 
                                  z[t]],z[t]]+F[kx[t], 
                          lx[Pn[t], ϕ1[t], ϕ2[t], kx[t], z[t]], 
                          z[t]]/.If[t < ts, gminsol1, gminsol2]]},{t, -2, tplot},PlotRange->{1.05, 1.61},AxesOrigin -> {-2, 1.05},AxesLabel -> {"t", Y},PlotStyle -> {{Thickness[ .007], 
                  Dashing[{ .05, 0}]}},DisplayFunction -> Identity];
pb=Plot[{If[t < 0,(Subscript[b, o]),b[t]/.If[t < ts, gminsol1, gminsol2]]},{t, -2, tplot},PlotRange -> {- .4, 1. },AxesOrigin -> {-2, 0},AxesLabel -> {"t", b},PlotStyle -> {{Thickness[ .007], 
                  Dashing[{ .05, 0}]}},DisplayFunction -> Identity];Show[GraphicsArray[{{pz,pEx,pPn},{pm,pπ,pχ},{plx,pln,pQn},{pix,pkx,pQx},{pY,pb}}]]

Labels={"m","Pn","Ex","ix","b","kx","z","Y","Qx","Qn","lx","ln","π","χ"};InitialSteadyState={(Subscript[m, o]),(Subscript[Pn, o]),(Subscript[Ex, o]),δ k  (Subscript[x, o]),(Subscript[b, o]),(Subscript[kx, o]),(Subscript[z, o]),(Subscript[Y, o]),(Subscript[Qx, o]),(Subscript[Qn, o]),(Subscript[lx, o]),lbar - (Subscript[lx, o]),(Subscript[π, o]),(Subscript[π, o])}//N;FinalSteadyState={(Subscript[m, f]),(Subscript[Pn, f]),(Subscript[Ex, f]),δ k  (Subscript[x, f]),(Subscript[b, f]),(Subscript[kx, f]),(Subscript[z, f]),(Subscript[Qx, f]) + (Subscript[Pn, o]) (Subscript[Qn, f]),(Subscript[Qx, f]),(Subscript[Qn, f]),(Subscript[lx, f]),lbar - (Subscript[lx, f]),(Subscript[π, f]),(Subscript[π, f])}//N;
(Dat=Table[{t,m[t],Pn[t],Ex[Pn[t], ϕ1[t]],ix[Pn[t], ϕ1[t], ϕ2[t], 
                                kx[t]],b[t],kx[t],z[t],(Subscript[Pn, o])G[(Subscript[kn, o]),lbar - 
                                        lx[Pn[t], ϕ1[t], ϕ2[t], 
                                        kx[t], z0],z0]+F[kx[t], 
                                  lx[Pn[t], ϕ1[t], ϕ2[t], kx[t], 
                                    z0], z[t]],F[kx[t], 
                                lx[Pn[t], ϕ1[t], ϕ2[t], kx[t], z0], 
                                z[t]],G[(Subscript[kn, o]),lbar - 
                                    lx[Pn[t], ϕ1[t], ϕ2[t], kx[t], 
                                      z0],z0],lx[Pn[t], ϕ1[t], ϕ2[t], 
                                kx[t], z0],lbar - 
                                lx[Pn[t], ϕ1[t], ϕ2[t], kx[t], 
                                  z0],χ[m[t], t]+((((((Subscript[k, 1])))^β) 
                                        Pn[t]^(1 - β))/(((((Subscript[k, 0])))^β) + ((((Subscript[k, 1])))^β) 
                                        Pn[t]^(1 - β)))((Pndot /. endovarsubs/.{Pn->Pn[t],ϕ1->ϕ1[t],ϕ2->ϕ2[t],b->b[t],kx->kx[t],z->z[t]})/(Pn[
                                        t]))/.If[t < ts, gminsol1, gminsol2]/.If[t < ts, lminsol1, lminsol2],χ[m[t], t]}/.If[t < ts, lminsol1, lminsol2]/.If[t < ts, gminsol1, gminsol2],{t,0,T,1/4}])                       //MatrixForm;
(Dat = 
                  Join[{Join[{"Initial Steady State"}, 
                        InitialSteadyState], 
                      Join[{"Final Steady State"}, FinalSteadyState], 
                      Join[{Time}, Labels]}, Dat]) // N // MatrixForm;
line={t,m[t],Pn[t],Ex[Pn[t], ϕ1[t]],ix[Pn[t], ϕ1[t], ϕ2[t], kx[t]],b[t],kx[t],z[t],(Subscript[Pn, o])G[(Subscript[kn, o]),lbar - 
                                  lx[Pn[t], ϕ1[t], ϕ2[t], kx[t], 
                                    z0],z0]+F[kx[t], 
                            lx[Pn[t], ϕ1[t], ϕ2[t], kx[t], z0], 
                            z[t]],F[kx[t], 
                          lx[Pn[t], ϕ1[t], ϕ2[t], kx[t], z0], 
                          z[t]],G[(Subscript[kn, o]),lbar - 
                              lx[Pn[t], ϕ1[t], ϕ2[t], kx[t], z0],z[t]],lx[Pn[t], ϕ1[t], ϕ2[t], kx[t], 
                          z0],lbar - 
                          lx[Pn[t], ϕ1[t], ϕ2[t], kx[t], z0],χ[m[t], t]+((((((Subscript[k, 1])))^β) 
                                      Pn[t]^(1 - β))/(((((Subscript[k, 0])))^β) + ((((Subscript[k, 1])))^β) 
                                        Pn[t]^(1 - β)))((Pndot /. endovarsubs/.{Pn->Pn[t],ϕ1->ϕ1[t],ϕ2->ϕ2[t],b->b[t],kx->kx[t],z->z[t]})/(Pn[
                                    t]))/.If[t < ts, gminsol1, gminsol2]/.If[t < ts, lminsol1, lminsol2],χ[m[t], t]}/.lminsol1/.gminsol1/.t->ts-(10^(-5));
Join[Part[Dat,Table[i,{i,1,(ts/(1/4))+1+2}]],{line},Part[Dat,Table[i,{i,(ts/(1/4))+1+2+1,(T/(1/4))+1+2+1}]]]//N//MatrixForm