Example 3: 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. There is capital accumulation in both export and nontraded sectors. The economy experiences a gradual rise in productivity. In this case, like Example 2, one has to solve the production and the consumption blocks together. Finally, one needs to solve for m by forward shooting.
THE OUTPUT OF THIS FILE CORRESPONDS TO FIGURE 10 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, f]) = (6/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]) + (Subscript[kn, o])) ; 
      eq2 :=  (Subscript[Qn, o]) == 
          (Subscript[Cn, o]) + ((α  (Subscript[Pk, o]) δ ((Subscript[kx, o]) + (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[kn, 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

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]];
in[Pn_, ϕ1_, ϕ3_, kn_] = 
              in /.  Solve[ϕ3 == ϕ1 Pk[
                          Pn] (1 + ψ'[
                              in/kn]), {in}][[1]];
lx[Pn_, ϕ1_, ϕ2_, ϕ3_, kx_, kn_, z_] = 
              lx /.  Solve[
                    Cn[Pn, Ex[
                            Pn, ϕ1]] + 
((α Pk[Pn] (ix[Pn, ϕ1, ϕ2, kx] + ψ[ix[Pn, ϕ1, ϕ2, kx]/kx] kx + in[Pn, ϕ1, ϕ3, kn] + ψ[in[Pn, ϕ1, ϕ3, kn]/kn] kn))/Pn) == G[kn, lbar - lx, z], {lx}][[1]];Pn[ϕ1_?NumberQ, ϕ2_?NumberQ, ϕ3_?
                  NumberQ, kx_?NumberQ, kn_?NumberQ, z_?NumberQ]:=Block[{Pn},{Pn}={Pn}/.FindRoot[{Pn (Subscript[G, 2])[kn, 
                                    lbar - 
                                      lx[Pn, ϕ1, ϕ2, ϕ3, kx, 
                                        kn, z], z]==(Subscript[F, 2])[kx, 
                                  lx[Pn, ϕ1, ϕ2, ϕ3, kx, kn, 
                                    z], z]},{Pn, Pnguess},MaxIterations -> 1000,WorkingPrecision -> digits];Pnguess = Pn;Pn];
endovarsubs={Pk -> Pk[Pn], Ex -> Ex[Pn, ϕ1], 
                ix -> ix[Pn, ϕ1, ϕ2, kx], 
                in -> in[Pn, ϕ1, ϕ3, kn], 
                lx -> lx[Pn, ϕ1, ϕ2, ϕ3, kx, kn, z]};

Pndot=(Dt[Pn]/.Solve[{Dt[ϕ2 == ϕ1 Pk[
                                      Pn] (1 + ψ'[ix/kx])],Dt[ϕ3 == ϕ1 Pk[
                                      Pn] (1 + ψ'[in/kn])],Dt[Pn (Subscript[G, 2])[kn, lbar - lx, z]==(Subscript[F, 2])[kx, lx, z]],
Dt[
                                Cn[Pn, 
                                      Ex[Pn, ϕ1]] + ((α Pk[Pn] (ix + ψ[ix/kx] kx + in + ψ[in/kn] kn))/Pn) == 
                                  G[kn, lbar - lx, z]]},{Dt[Pn], Dt[lx], Dt[ix], Dt[in]}][[1]])/.{Dt[ϕ1] -> ϕ1dot, 
                  Dt[ϕ2] -> ϕ2dot, 
                  Dt[ϕ3] -> ϕ3dot, Dt[kx] -> kxdot, 
                  Dt[kn] -> kndot, Dt[z] -> zdot};
ϕ1dot  = 
            0; ϕ2dot =  (ρ + δ) ϕ2 - ϕ1 ((Subscript[F, 1])[kx, lx, z] + 
                    Pk (ψ'[ix/kx] ix/kx - ψ[
                            ix/kx])); ϕ3dot =  (ρ + δ) ϕ3 - ϕ1 ((Subscript[G, 1])[kn, lbar - lx, z] + 
                    Pk (ψ'[in/kn] in/kn - ψ[in/kn])); 
          bdot =  Pn G[kn, lbar - lx, z] + F[kx, lx, z] + ρ b - Ex - 
              Pk (ix + ψ[ix/kx] kx + in + ψ[in/kn] kn);
kxdot = ix - δ kx; kndot = in - δ kn; 
          zdot  =   s ((Subscript[z, f]) - z);
system={Pndot,ϕ1dot,ϕ2dot,ϕ3dot,bdot,kxdot,kndot,zdot}/.endovarsubs;

Pndot /. endovarsubs/.{ϕ1->(Subscript[V, 2])[(Subscript[Pn, o]),(Subscript[Ex, o])],Pn->(Subscript[Pn, o]),ϕ2->(Subscript[V, 2])[(Subscript[Pn, o]),(Subscript[Ex, o])]Pk[(Subscript[Pn, o])],ϕ3->(Subscript[V, 2])[(Subscript[Pn, o]),(Subscript[Ex, o])]Pk[(Subscript[Pn, o])],b->(Subscript[b, o]),kx->(Subscript[kx, o]),kn->(Subscript[kn, o]),z->(Subscript[z, o])}//N

(* 
          Set up system and initial conditions for the ODE solver *)  
sys:=Thread[{Pn'[t],ϕ1'[t],ϕ2'[t],ϕ3'[t],b'[t],kx'[t],kn'[t],z'[t]}== (system/.{Pn->Pn[t],ϕ1->ϕ1[t],ϕ2->ϕ2[t],ϕ3->ϕ3[t],b->b[t],kx->kx[t],kn->kn[t],z->z[t]}),List];
ini:=Thread[{Pn[0],ϕ1[0],ϕ2[0],ϕ3[0],b[0],kx[0],kn[0],z[0]}=={Pnin,ϕ1in,ϕ2in,ϕ3in,bin,kxin,knin,zin} ,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)|| in[Pn[t], ϕ1[t], ϕ3[t], kn[t]]>((δ (Subscript[kn, o]))/5)||in[Pn[t], ϕ1[t], ϕ3[t], kn[t]]<((2 δ (Subscript[kn, o]))/1)|| lx[Pn[t], ϕ1[t], ϕ2[t], ϕ3[t], 
                    kx[t], kn[t], z[t]]>(1/5)||lx[Pn[t], ϕ1[t], ϕ2[t], ϕ3[t], 
                    kx[t], kn[t], z[t]]<(8/10)||b[t]>-10 (Subscript[Y, o])||b[t]<10(Subscript[Y, o])|| kx[t]>((Subscript[kx, o])/5)||kx[t]<(4/2) (Subscript[kx, o])|| kn[t]>((Subscript[kn, o])/5)||kn[t]<(4/2)(Subscript[kn, o]); 
  
sysdiff[t1_, t2_]:={Max[Abs[Table[(Ex[Pn[t], ϕ1[t]]/.gminsol)-(Ex[Pn[t], ϕ1[t]]/.gminsol1),{t, t1, t2, 1/10}]]],Max[Abs[Table[(ix[Pn[t], ϕ1[t], ϕ2[t], 
                                  kx[t]]/.gminsol)-(ix[Pn[t], ϕ1[t], ϕ2[t], 
                                  kx[t]]/.gminsol1),{t, t1, t2, 1/10}]]],Max[Abs[Table[(in[Pn[t], ϕ1[t], ϕ3[t], 
                                  kn[t]]/.gminsol)-(in[Pn[t], ϕ1[t], ϕ3[t], 
                                  kn[t]]/.gminsol1),{t, t1, t2, 1/10}]]],
Max[Abs[Table[(b[t]/.gminsol)-(b[t]/.gminsol1),{t, t1, t2, 1/10}]]],Max[Abs[Table[(kx[t]/.gminsol)-(kx[t]/.gminsol1),{t, t1, t2, 1/10}]]],Max[Abs[Table[(kn[t]/.gminsol)-(kn[t]/.gminsol1),{t, t1, t2, 1/10}]]],Max[Abs[Table[(z[t]/.gminsol)-(z[t]/.gminsol1),{t, t1, t2, 1/10}]]]};

distance[sol_, tsol_]:=
Block[{dis, dislast, tdis, steps, tbeg, tend, intsz},dis =  (10^7);dislast = (10^8);steps = Max[20, 4 Floor[tsol]] ;tbeg = 0 5  tsol/10;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]),in[
                                        Pn[t], ϕ1[t], ϕ3[t], 
                                        kn[t]]-δ (Subscript[kn, f]),b[t]-(Subscript[b, f]),kx[t]-(Subscript[kx, f]),kn[t]-(Subscript[kn, f]),z[t]-(Subscript[z, 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), tbeg];tend = Min[tdis + (intsz/steps), tsol];intsz = tend - tbeg;steps = 40];{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!";

SPHERE SEARCH

(* 
          Set up for solving for the global nonlinear saddle path *) 
Clear[compdis];
compdis[p_, epsilon_, theta_, phibeg_]:=Block[{Pnin,ϕ1in,ϕ2in,ϕ3in,bin,kxin,knin,zin,Pnmin1,ϕ1min1,ϕ2min1,ϕ3min1,ϕ1min2,Pnmin2,ϕ2min2,ϕ3min2,mindis,j1,j2,sol,tsol,dis,tdis,tmin,jmin,minsol,tmin2,jmin2,minsol2,mindis2},{ntraj, mindis, mindis2} = {ntrajini, (10^7), 
                      (10^7)};{(Subscript[kx, f]), (Subscript[kn, f]), (Subscript[lx, f])} = {(Subscript[kx, o]), (Subscript[kn, o]), (Subscript[lx, o])};
solve:=(phi=phibeg+((2π j)/ntraj);{ϕ1in,ϕ2in,ϕ3in}=p+epsilon{1Cos[
                                      theta]Cos[phi],1Sin[
                                      theta]Cos[phi],1Sin[
                                      phi]};{bin,kxin,knin,zin}={(Subscript[b, o]),(Subscript[kx, o]),(Subscript[kn, o]),(Subscript[z, o])};
   
Pnguess = (Subscript[Pn, o]);Pnin=Pn[ϕ1in, ϕ2in, ϕ3in, 
                              kxin, knin, zin];
  
ϕ(Subscript[1, f]) = ϕ1in;{kxguess, knguess, lxguess} = {(Subscript[kx, f]), 
                              (Subscript[kn, f]), (Subscript[lx, f])};{(Subscript[kx, f]), (Subscript[kn, f]), (Subscript[lx, f])} =. ;(Subscript[Pn, f]) = 
                            (Subscript[F, 2])[(Subscript[kx, f]), (Subscript[lx, f]), (Subscript[z, f])]/(Subscript[G, 2])[(Subscript[kn, f]), lbar - (Subscript[lx, f]), (Subscript[z, f])];
{(Subscript[kx, f]), (Subscript[kn, f]), (Subscript[lx, f])}={(Subscript[kx, f]), (Subscript[kn, f]), (Subscript[lx, f])}/.FindRoot[{Cn[(Subscript[Pn, f]), 
                                        Ex[(Subscript[Pn, f]), ϕ(Subscript[1, f])]]+α Pk[(Subscript[Pn, f])] δ ((Subscript[kx, f])+(Subscript[kn, f]))/(Subscript[Pn, f])==G[(Subscript[kn, f]), lbar - (Subscript[lx, f]), 
                                        (Subscript[z, f])],(ρ + δ) 
                                        Pk[(Subscript[Pn, f])] == 
                                        (Subscript[F, 1])[(Subscript[kx, f]), (Subscript[lx, f]), (Subscript[z, f])],(ρ + δ) 
                                        Pk[(Subscript[Pn, f])] == 
                                        (Subscript[G, 1])[(Subscript[kn, f]), lbar - (Subscript[lx, f]), (Subscript[z, f])]},{(Subscript[kx, f]), kxguess},{(Subscript[kn, f]), knguess},{(Subscript[lx, f]), lxguess},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[3, f]) = ϕ(Subscript[1, f]) Pk[(Subscript[Pn, f])];;(Subscript[b, f]) = (((Subscript[Ex, f]) + Pk[(Subscript[Pn, f])] δ ( (Subscript[kx, f]) + (Subscript[kn, f])) - F[(Subscript[kx, f]), (Subscript[lx, f]), (Subscript[z, f])] - (Subscript[Pn, f]) G[(Subscript[kn, f]), lbar - (Subscript[lx, f]), (Subscript[z, f])])/ρ);{Pn[0],ϕ1[0],ϕ2[0],ϕ3[0],b[0],kx[0],kn[0],z[0]}={Pnin,ϕ1in,ϕ2in,ϕ3in,bin,kxin,knin,zin} ;If[(box /. t -> 0) ≠  
                            True, Print["Out of Box."]; {dis, 
                              tdis} = {(10^7), 100}; Goto[skip]];{Pn[0],ϕ1[0],ϕ2[0],ϕ3[0],b[0],kx[0],kn[0],z[0]}=.;sol= First[NDSolve[{sys, ini},{Pn, ϕ1, ϕ2, ϕ3, 
                                    b, kx, kn, z},{t, 0, Tmax},WorkingPrecision -> digits ,MaxSteps -> 300,\(Method -> {"EventLocator", \ "Event" -> If[box == True, True, 
                                        False], "EventAction" :> \ 
                                      Throw[
                                        Null, \ "StopIntegration"], "EventLocationMethod" -> "StepBegin"}\)]];
tsol = 
                          sol[[1, 2, 1, 1, 2]];{dis, tdis} = distance[sol, tsol];Label[skip];(*Print[{j,dis,tdis,sol[[1, 2, 1, 1, 2]],ntraj,{ϕ1in,ϕ2in,ϕ3in}//N}];*));

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;dissol[j] = sol;jp[j]={ϕ1in,ϕ2in,ϕ3in},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]];Print[""];
\(For[ii = 1, 
                    ii ≤ Length[minpos], \(ii++\), 
                    Print["ii = ", 
                      ii, "   {mindis, tdis, tsol, jmin, ntraj} = ", {jdis[minpos[\([ii]\)]], 
                        jtdis[minpos[\([ii]\)]] // 
                          N, \((dissol[minpos[\([ii]\)]])\)[\([1, 2, 1, 1, 
                            2]\)], minpos[\([ii]\)], ntraj}]]\);If[
                    Length[minpos] == 2  &&  iflag == 1, 
                    trajtrack = 1];firstiter = 0;

For[ii = 1,ii ≤ Length[minpos],ii++,{mindis, tmin, jmin, ntraj, 
                            minsol} = {jdis[minpos[[ii]]], 
                            jtdis[minpos[[ii]]], minpos[[ii]], 
                            ntrajini, dissol[minpos[[ii]]]};{ϕ1min1,ϕ2min1,ϕ3min1}=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;minsol = sol;{ϕ1min1,ϕ2min1,ϕ3min1}={ϕ1in,ϕ2in,ϕ3in};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;minsol2 = minsol;{ϕ1min2,ϕ2min2,ϕ3min2}={ϕ1min1,ϕ2min1,ϕ3min1}];
  ];
  Print["Ring = ",jb,"    Minimum Distance = ",mindis2 // N,"    Minimum point = ",SetPrecision[{ ϕ1min2, ϕ2min2, ϕ3min2},precision],"    No. of Subdivisions of ϵ-ring > ",jmin2];Label[end1];{{mindis2, tmin2, minsol2},{ϕ1min2,ϕ2min2,ϕ3min2}}];

Clear[compdis2];
compdis2[p_, epsilon_, thetabeg_, phibeg_]:=Block[{ϕ1bmin,ϕ2bmin,ϕ3bmin,ϕ1bmin1,ϕ2bmin1,ϕ3bmin1,ϕ1bmin2,ϕ2bmin2,ϕ3bmin2,bmindis,jb1,jb2,sol,dis,tdis,btmin,jbmin,bminsol,btmin2,jbmin2,bminsol2,bmindis2,theta,jbdis,jbtdis,jbdissol,jbp},{nrings, bmindis, bmindis2} = {nringsini, (10^7), 
                  (10^7)};Label[begin];

For[jb = 1,jb ≤  nrings ,jb++,theta= thetabeg+((1π jb)/nrings);{{dis, tdis, sol},{ϕ1bmin,ϕ2bmin,ϕ3bmin}}=compdis[p, epsilon, theta, phibeg];
jbdis[jb] = dis;jbtdis[jb] = tdis;jbdissol[jb] = sol;jbp[jb]={ϕ1bmin,ϕ2bmin,ϕ3bmin}];

btab = 
                Table[jbdis[Mod[ib - 1, nrings] + 1] - 
                    jbdis[Mod[ib, nrings] + 1], {ib, 1, nrings}];bminpos = 
                Flatten[
                  Position[
                    Table[btab[[Mod[ib - 2, nrings] + 1]] > 0  && 
                        btab[[Mod[ib - 1, nrings] + 1]] < 0, {ib, 1, 
                        nrings}], True]];Print[" "];
For[ib = 1, 
                ib ≤ Length[bminpos], ib++, 
                Print["Initial Ring = ", 
                  ib, "   {bmindis, btdis, bjmin, nrings} = ", {jbdis[
                      bminpos[[ib]]], jbtdis[bminpos[[ib]]] // N, 
                    bminpos[[ib]], nrings}]];Print[" "];

For[ib = 1,ib ≤ Length[bminpos],ib++,{bmindis, btmin, jbmin, nrings, 
                        bminsol} = {jbdis[bminpos[[ib]]], 
                        jbtdis[bminpos[[ib]]], bminpos[[ib]], 
                        nringsini, jbdissol[bminpos[[ib]]]};
{jb1, jb2} = 
                      2 {jbmin - 1, jbmin + 1};{nrings, jbmin} = 2 {nrings, jbmin};{ ϕ1bmin1,  ϕ2bmin1,  ϕ3bmin1}=jbp[bminpos[[ib]]];Label[begin];

For[jb = jb1 + 1,jb ≤  jb2  ,jb += 2,theta= thetabeg+((1π jb)/nrings);{{dis, tdis, sol},{ ϕ1bmin, ϕ2bmin, ϕ3bmin}}=compdis[p, epsilon, theta, phibeg];
If[ Abs[bmindis - dis] ≤ 
                                Max[(tol/100), (gmindis/100)], Goto[end]; ];
If[dis ≤ bmindis,bmindis = dis;btmin = tdis;jbmin = jb;bminsol = sol;{ ϕ1bmin1,  ϕ2bmin1,  ϕ3bmin1}={ ϕ1bmin,  ϕ2bmin,  ϕ3bmin};Break[];];];

{jb1, jb2} = 
                      2 {jbmin - 1, jbmin + 1};{nrings, jbmin} = 2 {nrings, jbmin}; Goto[begin];Label[end];
If[bmindis < bmindis2,bmindis2 = bmindis;btmin2 = btmin;jbmin2 = jbmin;bminsol2 = bminsol;{ ϕ1bmin2,  ϕ2bmin2,  ϕ3bmin2}={ ϕ1bmin1, ϕ2bmin1, ϕ3bmin1}];
  ];Print[""];

Print["Ball = ",bit,"    Minimum Distance = ",bmindis2 // N,"    Minimum point = ",SetPrecision[{ ϕ1bmin2, ϕ2bmin2, ϕ3bmin2},precision],"    No. of Subdivisions of ϵ-ball > ",jbmin2];{{bmindis2, btmin2, bminsol2},{ ϕ1bmin2, ϕ2bmin2, ϕ3bmin2}}];

compdis3[p_, ϵ_]:=Block[{ϕ1gmin,ϕ2gmin,ϕ3gmin,dis,tdis,sol},If[
                firstsim == 
                  1, {pmin, epsilon} = {p, ϵ}, {pmin, 
                    epsilon} = (<< solution2)[[{3, 4}]]; iflag = 1];gmindis =  (10^7);iflag = 0;bit = 1;rinfo = {};Label[begin];Print[" "];Print["BEGINNING BALL ITERATION ",bit," WITH CENTER OF ϵ-ball at ",SetPrecision[pmin, precision]];{{dis, tdis, sol},{ϕ1gmin,ϕ2gmin,ϕ3gmin}}=compdis2[pmin, 
                  epsilon, -(π/1), -(π/2)];
If[dis ≤ gmindis,gmindis = dis;gtmin = tdis;gminsol = sol;pmin={ϕ1gmin,ϕ2gmin,ϕ3gmin};epsilon=(epsilon/1),epsilon = (epsilon/2)];
If[
                gtmin > T + 5  &&  iflag == 0, 
                iflag = 1; {gminsol, gtmin, pmin, 5  epsilon} >> solution1];
If[gtmin > 
                      Max[T + 5 , (<< solution1)[[2]] + 5] &&  
                    iflag == 1 ,{gminsol1, 
                        gtmin1} = (<< solution1)[[{1, 2}]];If[Max[sysdiff] > tol,Print[
                          sysdiff]; {gminsol, gtmin, pmin, 5  epsilon} >> 
                          solution1,{gminsol, gtmin, pmin, 5  epsilon} >> 
                            solution1;Print[sysdiff];Goto[end]];];
Print["AT THE END OF BALL ITERATION = ",bit,":    CURRENT MINIMUM DISTANCE = ",gmindis // N,"     TIME OF CLOSEST APPROACH = ",gtmin];
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];];bit++;
If[gmindis < tol1,Print["Target Tolerance = ", 
                      tol, "   Tolerance Achieved = ", gmindis]; Goto[end]; ];
Goto[begin];Label[end];{gtmin, gminsol}];

(* INITIALIZATION CELL - Initialization for finding the required trajectory. *)
{(Subscript[ϕ1, o]),(Subscript[ϕ2, o]),(Subscript[ϕ3, o])}={(Subscript[V, 2])[(Subscript[Pn, o]), (Subscript[Ex, o])], (Subscript[V, 2])[(Subscript[Pn, o]), (Subscript[Ex, o])] Pk[(Subscript[Pn, o])], 
                (Subscript[V, 2])[(Subscript[Pn, o]), (Subscript[Ex, o])] Pk[(Subscript[Pn, o])]};          
T = 20;                   
tol=(10^(-4));           ;
epsilon= (((10^(-0)))/400)(((((Subscript[ϕ1, o])^2)+((Subscript[ϕ2, o])^2)+((Subscript[ϕ3, o])^2)))^((1/2)));                                                  ;
pini = {33.7199760705672728487839966976699539261890430035449351195041420341421690805500453662782209675928149344849629916688869`99.42628808407845, 
              35.4249584668686313371936988183827929865178420087345698083919776253922422506501228122305480511066973910650226068006024`99.35048033187793, 
              37.4893935449129161365845883172308259875665710253640411224433310611486987169727526542323362851396662530268295772234584`99.38367224501347};

Timing[precision=16;firstsim = 1;Tmax = 200;ϵmax = epsilon;pinit = pini;gmindis = (10^7);trajtrack = 0;
While[True,localmin = 0;Label[begin];$MinPrecision = digits = precision;ntrajini = 20;nringsini = 7;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, 0, 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} = 
                compdis3[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,
{gminsol1, 
                        gtmin1} = (<< solution2)[[{1, 2}]];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[end1]]];];   ;
Label[end1];{gminsol1, gtmin1} = (<< solution1)[[{1, 2}]]; {gminsol, gtmin} = (<< solution2)[[{1, 2}]];
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."]]]

<< solution1

LINE SEARCH for m

{gminsol, gtmin, 
              pmin} = (<< solution1)[[{1, 2, 3}]];
ϕ(Subscript[1, f]) = pmin[[1]];(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, f]), 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[kn, f]), (Subscript[lx, f])} =. ;{(Subscript[kx, f]), (Subscript[kn, f]), (Subscript[lx, f])}={(Subscript[kx, f]), (Subscript[kn, f]), (Subscript[lx, f])}/.FindRoot[{Cn[(Subscript[Pn, f]), Ex[(Subscript[Pn, f]), ϕ(Subscript[1, f])]]+α Pk[(Subscript[Pn, f])] δ ((Subscript[kx, f])+(Subscript[kn, f]))/(Subscript[Pn, f])==G[(Subscript[kn, f]), lbar - (Subscript[lx, f]), (Subscript[z, f])],(ρ + δ) Pk[(Subscript[Pn, f])] == 
                            (Subscript[F, 1])[(Subscript[kx, f]), (Subscript[lx, f]), (Subscript[z, f])],(ρ + δ) Pk[(Subscript[Pn, f])] == 
                            (Subscript[G, 1])[(Subscript[kn, f]), lbar - (Subscript[lx, f]), (Subscript[z, f])]},{(Subscript[kx, f]),(Subscript[kx, o])},{(Subscript[kn, f]),(Subscript[kn, 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[3, f]) = ϕ(Subscript[1, f]) Pk[(Subscript[Pn, f])];; 
        (Subscript[b, f]) = (((Subscript[Ex, f]) + Pk[(Subscript[Pn, f])] δ ( (Subscript[kx, f]) + (Subscript[kn, f])) - F[(Subscript[kx, f]), (Subscript[lx, f]), (Subscript[z, f])] - (Subscript[Pn, f]) G[(Subscript[kn, f]), 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];
{{(Subscript[m, o]), (Subscript[Ex, o]), (Subscript[Pk, o]), (Subscript[Pn, o]), (Subscript[kn, 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])}, {(Subscript[m, f]), (Subscript[Ex, f]), (Subscript[Pk, f]), (Subscript[Pn, f]), (Subscript[kn, f]), 
                (Subscript[kx, f]), (Subscript[ln, f]), (Subscript[lx, f]), (Subscript[Y, f]), (Subscript[Qn, f]), (Subscript[Qx, f]), (Subscript[b, f]), g, 
                (Subscript[z, f]), (Subscript[π, f]), (Subscript[Ex, f])/(Subscript[m, f]), (Subscript[m, f])/(Subscript[Y, f])}} // N // 
          MatrixForm

Clear[pi];χ[m_?NumberQ, 
              t_?NumberQ] := (( (Subscript[ϕ, 2])[Pn[t], m])/( ϕ(Subscript[1, f]))) - ρ /. gminsol

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, gmindis] (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] /. lminsol) - (m[t] /. lminsol1), {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[0, 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];Print["\nBEGINNING ITERATION ",lit," WITH CENTER OF ϵ-segment at ",{pmin, epsilon}];
For[j = -1,j ≤  1  ,j = j + 2,sol = First[NDSolve[{m'[t] == 
                                        g - χ[m[t], t] m[t] ,m[0] == 
                                        pmin + j epsilon},{m},{t, 0, gtmin},\(Method -> {"EventLocator", \ "Event" -> m[t] > Min[m\_o, m\_f]\/2 || 
                                        m[t] < 
                                        2  Max[m\_o, 
                                        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] = pmin + j epsilon;jm[j] = m[sol[[1, 2, 1, 1, 2]]] /. sol;];
If[
                    jm[1] > (Subscript[m, f])  &&  
                      jm[-1] > (Subscript[m, f]), {lmindis, ltmin, lminsol, 
                        pmin} = {jdis[-1], jtdis[-1], dissol[-1], 
                        jp[-1]}, 
                    If[jm[1] > (Subscript[m, f]), 
                      epsilon = 
                        epsilon/2, {lmindis, ltmin, lminsol, 
                          pmin} = {jdis[1], jtdis[1], dissol[1], jp[1]}]];
Print["AT THE END OF ITERATION = ",lit,"     TIME OF CLOSEST APPROACH = ",ltmin];If[lit == 300, Abort[]];

If[
                    ltmin > T + 10  &&  iflag == 0, 
                    iflag = 1; {lminsol, ltmin, lmindis, pmin} >> solutionm1; 
                    Print[""]; Print[step2];];
If[ltmin > 
                          Max[T + 10 , (<< solutionm1)[[2]] + 
                              5] &&  iflag == 1 ,{lminsol1, 
                            ltmin1} = (<< solutionm1)[[{1, 2}]];If[sysdiffm[0, T] > tol,Print[""]; difftminm; 
                            Print[notsatistminm]; {lminsol, ltmin, lmindis, 
                                pmin} >> solutionm1,Print[""];difftminm;Print[satistminm];{lminsol, ltmin, lmindis, pmin} >> 
                                solutionm1;Goto[end]]];lit++;Goto[begin];Label[end];

  ;
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."]]
{lminsol, ltmin, lmindis, 
                      pmin}];

Timing[precision=16;Tmax = gtmin;$MinPrecision = digits = precision;
sol=First[NDSolve[{m'[t] == g - χ[m[t], t] m[t] ,m[0] == (Subscript[m, o])},{m},{t, 0, Tmax},MaxSteps -> 1000,WorkingPrecision -> digits]];tsol = sol[[1, 2, 1, 1, 2]];lmindis = 1;{lmindis, ltmin} = distancem[sol, tsol];Print[{tsol, lmindis, ltmin}];
  ;
{lminsol, ltmin, lmindis, pmin}=compdism[(Subscript[m, o]),((Subscript[m, o])/4)];Print[congratulationsm]; ]

RESULTS

(*temp2=<<solution1*)

{temp1, temp2} = << fig12; T = 20

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

Clear[pi];χ[m_?NumberQ, 
              t_?NumberQ] := (( (Subscript[ϕ, 2])[Pn[t], m])/( ϕ(Subscript[1, f]))) - ρ /. gminsol

{lminsol, 
          gminsol} =   {temp1[[1]], 
          temp2[[1]]};
tplot  = T + 5;
pm=Plot[{(Subscript[m, o]),m[t]/.lminsol,(Subscript[m, f])},{t, 0, tplot},PlotRange -> All,AxesLabel->{"",m},DisplayFunction -> Identity];
pEx=Plot[{(Subscript[Ex, o]),Ex[Pn[t], ϕ1[t]]/.gminsol,(Subscript[Ex, f])},{t, 0, tplot},PlotRange -> All,AxesLabel->{"",Ex},DisplayFunction -> Identity];
pPn=Plot[{(Subscript[Pn, o]),Pn[t]/.gminsol,(Subscript[Pn, f])},{t, 0, tplot},PlotRange -> All,AxesLabel->{"",Pn},DisplayFunction -> Identity];
pix=Plot[{δ (Subscript[kx, o]),ix[Pn[t], ϕ1[t], ϕ2[t], kx[t]]/.gminsol,δ (Subscript[kx, f])},{t, 0, tplot},PlotRange -> All,AxesLabel->{"",ix},DisplayFunction -> Identity];
pin=Plot[{δ (Subscript[kn, o]),in[Pn[t], ϕ1[t], ϕ3[t], kn[t]]/.gminsol,δ (Subscript[kn, f])},{t, 0, tplot},PlotRange -> All,AxesLabel->{"",in},DisplayFunction -> Identity];
pb=Plot[{(Subscript[b, o]),b[t]/.gminsol,(Subscript[b, f])},{t, 0, tplot},PlotRange -> All,AxesLabel->{"",b},DisplayFunction -> Identity];
pkx=Plot[{(Subscript[kx, o]),kx[t]/.gminsol,(Subscript[kx, f])},{t, 0, tplot},PlotRange -> All,AxesLabel->{"",kx},DisplayFunction -> Identity];
pkn=Plot[{(Subscript[kn, o]),kn[t]/.gminsol,(Subscript[kn, f])},{t, 0, tplot},PlotRange -> All,AxesLabel->{"",kn},DisplayFunction -> Identity];
pz=Plot[{(Subscript[z, o]),z[t]/.gminsol,(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], ϕ3[t], kx[t], 
                        kn[t], z[t]], z[t]]/.gminsol,(Subscript[Qx, f])},{t, 0, tplot},PlotRange -> All,AxesLabel->{"",Qx},DisplayFunction -> Identity];
pQn=Plot[{(Subscript[Qn, o]),G[kn[t], 
                      lbar - lx[Pn[t], ϕ1[t], ϕ2[t], ϕ3[t], 
                          kx[t], kn[t], z[t]], z[t]]/.gminsol,(Subscript[Qn, f])},{t, 0, tplot},PlotRange -> All,AxesLabel->{"",Qn},DisplayFunction -> Identity];
pY=Plot[{(Subscript[Y, o]),Pn[t]G[(Subscript[kn, o]),lbar - 
                              lx[Pn[t], ϕ1[t], ϕ2[t], ϕ3[t], 
                                kx[t], kn[t], z[t]],z[t]]+F[kx[t], 
                        lx[Pn[t], ϕ1[t], ϕ2[t], ϕ3[t], kx[t], 
                          kn[t], z[t]], z[t]]/.gminsol,(Subscript[Y, f])},{t, 0, tplot},PlotRange -> All,AxesLabel->{"",Y},DisplayFunction -> Identity];
plx=Plot[{(Subscript[lx, o]),lx[Pn[t], ϕ1[t], ϕ2[t], ϕ3[t], 
                      kx[t], kn[t], z[t]]/.gminsol,(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], ϕ3[t], kx[t], 
                        kn[t], z[t]]/.gminsol,lbar-(Subscript[lx, f])},{t, 0, tplot},PlotRange -> All,AxesLabel->{"",ln},DisplayFunction -> Identity];
pχ=Plot[{(Subscript[π, o]),χ[m[t], t]/.lminsol,(Subscript[π, f])},{t, 0, tplot},PlotRange -> All,AxesOrigin->{0,  .99  (Subscript[π, f])},AxesLabel->{"",π},DisplayFunction -> Identity];
Show[GraphicsArray[{{pm,pPn,pEx,pb},{pix,pin,pkx,pkn},{pz,pQx,pQn,pY},{plx,pln,pχ}}]]

tplot =T;
pz=Plot[{If[t < 0,(Subscript[z, o]),z[t]/.gminsol]},{t, -2, tplot},PlotRange -> { .95, 1.21},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]]))/.gminsol]},{t, -2, tplot},PlotRange->{ .245,  .295},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]/.gminsol]},{t, -2, tplot},PlotRange -> { .95, 1.25},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]/.lminsol]},{t, -2, tplot},PlotRange -> { .09,  .125},AxesOrigin -> {-2,  .09},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], ϕ3 -> ϕ3[t], 
                                  b -> b[t], kn -> kn[t], 
                                  kx -> kx[t], z -> z[t]}) /. 
                        gminsol /. lminsol]},{t, -2, tplot},PlotRange -> { .36,  .51},AxesOrigin->{-2,.36},AxesLabel -> {"t", π},PlotStyle -> {{Thickness[ .007], 
                  Dashing[{ .05, 0}]}},DisplayFunction -> Identity];pχ=Plot[{If[t < 0,(Subscript[π, o]),χ[m[t], t] /. gminsol /. lminsol]},{t, -2, tplot},PlotRange -> { .36,  .51},AxesOrigin->{-2,.36},AxesLabel -> {"t", χ},PlotStyle -> {{Thickness[ .007], 
                  Dashing[{ .05, 0}]}},DisplayFunction -> Identity];
pin=Plot[{If[t < 0,δ (Subscript[kn, o]),in[Pn[t], ϕ1[t], ϕ3[t], kn[t]]/.gminsol]},{t, -2, tplot},PlotRange -> { .095,  .125},AxesOrigin -> {-2,  .095},AxesLabel->{"t",(Subscript["I", n])},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]]/.gminsol]},{t, -2, tplot},PlotRange -> {},AxesOrigin -> {-2,  .455},AxesLabel->{"t",(Subscript[L, n])},PlotStyle -> {{Thickness[ .007], 
                      Dashing[{ .05, 0}]}},DisplayFunction -> Identity];*)
pkn=Plot[{If[t < 0,(Subscript[kn, o]),kn[t]/.gminsol]},{t, -2, tplot},PlotRange -> {1.98, 2.22},AxesOrigin -> {-2, 1.98},AxesLabel->{"t",(Subscript[K, 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], ϕ3[t], 
                              kn[t], kx[t], z[t]],z[t]]/.gminsol]},{t, -2, tplot},PlotRange -> { .595,  .685},AxesOrigin -> {-2,  .595},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]]/.gminsol]},{t, -2, tplot},PlotRange -> { .085,  .155},AxesOrigin -> {-2,  .085},AxesLabel->{"t",(Subscript["I", x])},PlotStyle -> {{Thickness[ .007], 
                  Dashing[{ .05, 0}]}},DisplayFunction -> Identity];
pkx=Plot[{If[t < 0,(Subscript[kx, o]),kx[t]/.gminsol]},{t, -2, tplot},PlotRange -> {2.05, 2.6},AxesOrigin -> {-2, 2.05},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], ϕ3[t], kn[t], 
                          kx[t], z[t]], z[t]]/.gminsol]},{t, -2, tplot},PlotRange -> { .6,  .92},AxesOrigin -> {-2,  .6},AxesLabel->{"t",(Subscript[Q, x])},PlotStyle -> {{Thickness[ .007], 
                  Dashing[{ .05, 0}]}},DisplayFunction -> Identity];plx=Plot[{If[t < 0,(Subscript[lx, o]),lx[Pn[t], ϕ1[t], ϕ2[t], ϕ3[t], 
                        kn[t], kx[t], z[t]]/.gminsol]},{t, -2, tplot},PlotRange -> { .48,  .6},AxesOrigin -> {-2,  .48},AxesLabel->{"t",(Subscript[L, 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], ϕ3[t], 
                                  kn[t], kx[t], z[t]],z[t]]+F[kx[t], 
                          lx[Pn[t], ϕ1[t], ϕ2[t], ϕ3[t], 
                            kn[t], kx[t], z[t]], z[t]]/.gminsol]},{t, -2, tplot},PlotRange -> {1.2, 1.6},AxesOrigin->{-2,1.2},AxesLabel -> {"t", Y},PlotStyle -> {{Thickness[ .007], 
                  Dashing[{ .05, 0}]}},DisplayFunction -> Identity];
pb=Plot[{If[t < 0,(Subscript[b, o]),b[t]/.gminsol]},{t, -2, tplot},PlotRange -> {-1.8, 0},AxesOrigin -> {-2, 0},AxesLabel -> {"t", b},PlotStyle -> {{Thickness[ .007], 
                  Dashing[{ .05, 0}]}},DisplayFunction -> Identity];
Show[GraphicsArray[{{pz,pEx,pPn},{pm,pπ,pχ},{pin,pkn,pQn},{pix,pkx,pQx},{plx,pY,pb}}]]

Labels={"m","Pn","Ex","ix","in","b","kx","kn","z","Y","Qx","Qn","lx","ln","pi"};InitialSteadyState={(Subscript[m, o]),(Subscript[Pn, o]),(Subscript[Ex, o]),δ (Subscript[kx, o]),δ (Subscript[kn, o]),(Subscript[b, o]),(Subscript[kx, o]),(Subscript[kn, o]),(Subscript[z, o]),(Subscript[Y, o]),(Subscript[Qx, o]),(Subscript[Qn, o]),(Subscript[lx, o]),lbar - (Subscript[lx, o]),(Subscript[π, o])}//N;FinalSteadyState={(Subscript[m, f]),(Subscript[Pn, f]),(Subscript[Ex, f]),δ (Subscript[kx, f]),δ (Subscript[kn, f]),(Subscript[b, f]),(Subscript[kx, f]),(Subscript[kn, f]),(Subscript[z, f]),(Subscript[Y, f]),(Subscript[Qx, f]),(Subscript[Qn, f]),(Subscript[lx, f]),lbar - (Subscript[lx, 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]],in[Pn[t], ϕ1[t], ϕ3[t], 
                                kn[t]],b[t],kx[t],kn[t],z[t],Pn[t] 
                                  G[kn[t], 
                                    lbar - 
                                      lx[Pn[t], ϕ1[t], ϕ2[
                                        t], ϕ3[t], kx[t], kn[t], z[t]], 
                                    z[t]]+F[kx[t], 
                                  lx[Pn[t], ϕ1[t], ϕ2[t], ϕ3[
                                      t], kx[t], kn[t], z[t]], z[t]],F[kx[t], 
                                lx[Pn[t], ϕ1[t], ϕ2[t], ϕ3[t], 
                                  kx[t], kn[t], z[t]], z[t]],G[kn[t], 
                                lbar - 
                                  lx[Pn[t], ϕ1[t], ϕ2[t], ϕ3[
                                      t], kx[t], kn[t], z[t]], z[t]],lx[
                                Pn[t], ϕ1[t], ϕ2[t], ϕ3[t], 
                                kx[t], kn[t], z[t]],lbar - 
                                lx[Pn[t], ϕ1[t], ϕ2[t], ϕ3[t], 
                                  kx[t], kn[t], z[t]],pi[m[t], t]}/.lminsol/.gminsol,{t,0,20,1/4}])//MatrixForm;
(Dat = 
                Join[{Join[{"Initial Steady State"}, InitialSteadyState], 
                    Join[{"Final Steady State"}, FinalSteadyState], 
                    Join[{Time}, Labels]}, Dat]) // N // 
          MatrixForm