This file implements the algorithm 'Forward Shoot 3D'  outlined in paper titled, "Smart Forward Shooting" which solves models with three jump variables by forward shooting.

The algorithm and the program written to implement the algorithm is free for use as long as proper credit is given to the authors of the paper and this first cell of the program is kept unchanged. This program is licensed under the CC-GNU GPL license with no expressed or implied guarantees.

Works with Mathematica Version 5.2 or later.

Version: April 07, 2008.
Author: Manoj Atolia

For details of the particular model being solved here, see Example 3 of the paper.

Off[General::spell];Off[General::spell1];Off[Unset::norep];Off[Unset::wrsym];Off[InterpolatingFunction::dmval];Off[Precision::precsm];

(* The precision of preliminary calculations. *)
{precision = 100, $MinPrecision=digits=precision};

Initialization

(* Set the directory where Mathematica will store intermediate results. *)
SetDirectory["c:"];

(* ALWAYS ENTER NUMERICAL VALUES AS FRACTIONS AND NOT AS DECIMALS TO ENSURE REQUISITE PRECISION IN CALCULATIONS *)   
{τ=(50/100),σ=(3/4),ρ=(1/10),β=(1/2),δ=(1/20),Ω=2}//N                                   (* Parameters of the model *)

{(Subscript[bbyY, o])=(25/100),(Subscript[mbyY, o])=(8/100),(Subscript[FbyY, o])=(8/100),(Subscript[π, o])=(25/100),(Subscript[θ, o])=(1/2)}//N        (* Data for calibration *)

{(Subscript[k, o])=1}//N                                                                          (* Normalization *)

{(Subscript[v, f])=(15/100),(Subscript[v, o])=(20/100),(Subscript[ω, 1])=(1/4),(Subscript[ω, 2])=(1/2)}//N                                                                 (* Shock - Exogenous rise in productivity *)

Steady State

(* Define utility function, production function etc. *)
U[c_]:=(1/(1-(1/τ)))(c^(1-(1/τ)));(Subscript[U, 1])[c_]:=Evaluate[D[U[c],c]];
ϕ[m_,F_]:=((Subscript[g, 0])/(1-(1/τ)))(((((((Subscript[g, 1])(m^(1-(1/σ)))+(1-(Subscript[g, 1]))(F^(1-(1/σ)))))^(σ/(σ-1)))))^(1-(1/τ)));(Subscript[ϕ, 1])[m_,F_]:=Evaluate[D[ϕ[m,F],m]];(Subscript[ϕ, 2])[m_,F_]:=Evaluate[D[ϕ[m,F],F]];
ψ[x_]:=(n/2)(((x-δ))^2);n=(1/(Ω δ));
Y[k_]:=((((Subscript[a, 0])+(Subscript[a, 1])(k^((β-1)/β))))^(β/(β-1)));(Subscript[Y, 1])[k_]:=Evaluate[D[Y[k],k]];
(* Avoid using symbols for subscripts for representing derivatives. Thus, (Subscript[Y, 1])[k,z] is (Subscript[Y, K])[k,z] of the paper.*)
(* Since we are using symbols o and f as subscripts for initial and final steady states avoid using o and f as names of variables or function. *)

(* Calibrate and compute the original steady state *)
(Subscript[Y, o])=Y[(Subscript[k, o])];{(Subscript[c, o])=(Subscript[Y, o])-δ (Subscript[k, o]),(Subscript[m, o])= (Subscript[mbyY, o]) (Subscript[Y, o]),(Subscript[F, o])= (Subscript[FbyY, o]) (Subscript[Y, o]),b=(Subscript[bbyY, o]) (Subscript[Y, o]),(Subscript[z, o])=(Subscript[π, o])(Subscript[m, o])+(Subscript[v, o])Y[(Subscript[k, o])]-ρ b};
{(Subscript[a, 0]),(Subscript[a, 1])}=.;{(Subscript[a, 0]),(Subscript[a, 1])}={(Subscript[a, 0]),(Subscript[a, 1])}/.FindRoot[{(Subscript[Y, 1])[(Subscript[k, o])] == ((ρ+δ)/(1-(Subscript[v, o]))),(Subscript[θ, o])==(((Subscript[Y, 1])[(Subscript[k, o])](Subscript[k, o]))/(Subscript[Y, o]))},{(Subscript[a, 0]),1},{(Subscript[a, 1]),1},MaxIterations->1000,AccuracyGoal->digits,WorkingPrecision->digits];{(Subscript[g, 0]),(Subscript[g, 1])}=.;{(Subscript[g, 0]),(Subscript[g, 1])}={(Subscript[g, 0]),(Subscript[g, 1])}/.FindRoot[{(((Subscript[ϕ, 1])[(Subscript[m, o]),(Subscript[F, o])])/((Subscript[U, 1])[(Subscript[c, o])])) == ρ+(Subscript[π, o]),(((Subscript[ϕ, 2])[(Subscript[m, o]),(Subscript[F, o])])/((Subscript[U, 1])[(Subscript[c, o])])) == ρ},{(Subscript[g, 0]),1},{(Subscript[g, 1]),(1/2)},MaxIterations->1000,AccuracyGoal->digits,WorkingPrecision->digits];

(* Compute the new steady state *)
(Subscript[k, f])=.;(Subscript[k, f])=(Subscript[k, f])/.FindRoot[{(Subscript[Y, 1])[(Subscript[k, f])]==((ρ+δ)/(1-(Subscript[v, f])))},{(Subscript[k, f]),(Subscript[k, o])},MaxIterations->1000,AccuracyGoal->digits,WorkingPrecision->digits];(Subscript[Y, f])=Y[(Subscript[k, f])];(Subscript[c, f])=(Subscript[Y, f])-δ (Subscript[k, f]);(Subscript[z, f])=(Subscript[z, o])+(Subscript[v, f])Y[(Subscript[k, f])]-(Subscript[v, o])Y[(Subscript[k, o])] ;
{(Subscript[m, f]),(Subscript[F, f]),(Subscript[π, f])}=.;{(Subscript[m, f]),(Subscript[F, f]),(Subscript[π, f])}={(Subscript[m, f]),(Subscript[F, f]),(Subscript[π, f])}/.FindRoot[{(((Subscript[ϕ, 1])[(Subscript[m, f]),(Subscript[F, f])])/((Subscript[U, 1])[(Subscript[c, f])])) == ρ+(Subscript[π, f]),(((Subscript[ϕ, 2])[(Subscript[m, f]),(Subscript[F, f])])/((Subscript[U, 1])[(Subscript[c, f])])) == ρ,(Subscript[z, f])==(Subscript[π, f])(Subscript[m, f])+(Subscript[v, f])Y[(Subscript[k, f])]-ρ b},{(Subscript[m, f]),(3/2)(Subscript[m, o])},{(Subscript[F, f]),(3/4)(Subscript[F, o])},{(Subscript[π, f]),(Subscript[π, o])},AccuracyGoal->digits,WorkingPrecision->digits,MaxIterations->1000];
{{"m","F","c","k","Y","b","z","pi","c/m","m/Y"},{(Subscript[m, o]),(Subscript[F, o]),(Subscript[c, o]),(Subscript[k, o]),(Subscript[Y, o]),b,(Subscript[z, o]),(Subscript[π, o]),(Subscript[c, o])/(Subscript[m, o]),(Subscript[m, o])/(Subscript[Y, o])},{(Subscript[m, f]),(Subscript[F, f]),(Subscript[c, f]),(Subscript[k, f]),(Subscript[Y, f]),b,(Subscript[z, f]),(Subscript[π, f]),(Subscript[c, f])/(Subscript[m, f]),(Subscript[m, f])/(Subscript[Y, f])}}//N//MatrixForm

Global Nonlinear Saddle Path

(* Solve for the endogenous variables as functions of the variables in the system. For consistency maintain same order of variables in arguments for functions of endogenous variables as in the system *)
{
rfunc[m_,c_,F_]:=(((Subscript[ϕ, 2])[m,F])/((Subscript[U, 1])[c]));
pifunc[m_,c_, F_]:=(((Subscript[ϕ, 1])[m,F]-(Subscript[ϕ, 2])[m,F])/((Subscript[U, 1])[c]));
};

(* Specify the functions, IF ANY, that, for a given value of t, returns the values of the (shock) variables that change exogenously. *)
z[t_]:=(Subscript[z, o])+((Subscript[z, f])-(Subscript[z, o]))(1-Exp[-(Subscript[ω, 1])t]);
v[t_]:=(Subscript[v, o])+((Subscript[v, f])-(Subscript[v, o]))(1-Exp[-(Subscript[ω, 2])t]);

(* Specify the nonlinear system with jump variables in the core dynamic system in the beginning and the state variables at the end *)
{
mdot = z+r b-v Y[k] - pi m; 
cdot=τ c (r-ρ);
λ1dot=(ρ-r)λ1;
λ2dot=(ρ+δ)λ2-((1-v)(Subscript[Y, 1])[k]+ψ'[(in/k)](in/k)-ψ[(in/k)])λ1;
indot = din/.Solve[Dt[λ2==(1+ψ'[(in/k)])λ1]/.{Dt[k]->kdot,Dt[in]->din,Dt[λ1]->λ1dot,Dt[λ2]->λ2dot},din][[1]]/.{λ1->(Subscript[U, 1])[c],λ2->(1+ψ'[in/k])(Subscript[U, 1])[c]};
kdot = in-δ k;
Fdot = Y[k] - c - in - (n/2)((((in/k)-δ))^2)k;
system  ={mdot,cdot,indot,kdot,Fdot};
};

(* INITIALIZATION CELL - Initialization for solving for the global saddle path. *)
{
  T = 20;                                                             (* solution time (years). BEWARE!!! Setting T larger than necessary can substantially increase computational time.*); 
tol = (10^(-4));                                                      (* desired accuracy of the solutions of the variables in the system *);
precision = 16;                                            (* starting value of precision of calculations *) ;
epsilon= ((10^(-0))/40)(((((Subscript[m, o])^2)+((Subscript[c, o])^2)))^(1/2));                    (* radius of the ϵ-ball *)
pini={(Subscript[m, o]),(Subscript[c, o]),δ (Subscript[k, o])};                                  (* center of the initial ϵ-ball *);
};
(* IN THE UNLIKELY EVENT THAT THE PROGRAM FAILS AT THE FIRST ATTEMPT TRY INCREASING THE VALUE OF 'precision'. ONE MAY ALSO EXPERIMENT WITH DIFFERENT VALUES OF 'epsilon' and 'pini'. *)

(* THE OUTPUT OF THIS CELL SHOWS THE SYNTAX FOR THE solve COMMAND *)
?solve

(* MAIN CELL - Solving for the global nonlinear saddlepath. *)
{gtmin,gminsol}=solve[system,{t},{m,c,in,k,F},{(Subscript[m, o]),(Subscript[c, o]),δ (Subscript[k, o]),(Subscript[k, o]),(Subscript[F, o])},{(Subscript[m, f]),(Subscript[c, f]),δ (Subscript[k, f]),(Subscript[k, f]),(Subscript[F, f])},{{r,rfunc,{m,c,F}},{pi,pifunc,{m,c,F}},{z,z,{t}},{v,v,{t}}},{},{},{T,tol,precision,epsilon},pini]
MaxMemoryUsed[]/1000000MB//N

(* Plot the final solutions. *)
tplot = T+10;
pc=Plot[{If[t<0,(Subscript[c, o]),c[t]]/.gminsol,(Subscript[c, f])},{t,-tplot/10,tplot},PlotRange->{.319,.334},AxesOrigin->{-tplot/10,.319},AxesLabel->{"Years",None},PlotStyle->{{Thickness[.007],Dashing[{.05,.0}]},{Thickness[.007],Dashing[{.02,0.02}]},{Thickness[.007],Dashing[{.01,0.01}]}},PlotLabel->"Real Consumption",DisplayFunction->Identity];
pin=Plot[{If[t<0,δ (Subscript[k, o]),in[t]]/.gminsol,δ (Subscript[k, f])},{t,-tplot/10,tplot},PlotRange->{.0495,.054},AxesOrigin->{-tplot/10,.0495},AxesLabel->{"Years",None},PlotStyle->{{Thickness[.007],Dashing[{.05,.0}]},{Thickness[.007],Dashing[{.02,0.02}]},{Thickness[.007],Dashing[{.01,0.01}]}},PlotLabel->"Investment",DisplayFunction->Identity];
pm=Plot[{If[t<0,(Subscript[m, o]),m[t]]/.gminsol,(Subscript[m, f])},{t,-tplot/10,tplot},PlotRange->{.018,.032},AxesOrigin->{-tplot/10,.018},AxesLabel->{"Years",None},PlotStyle->{{Thickness[.007],Dashing[{.05,.0}]},{Thickness[.007],Dashing[{.02,0.02}]},{Thickness[.007],Dashing[{.01,0.01}]}},PlotLabel->"Real Money Balances",DisplayFunction->Identity];
pF=Plot[{If[t<0,(Subscript[F, o]),F[t]]/.gminsol,(Subscript[F, f])},{t,-tplot/10,tplot},PlotRange->{.0295,.034},AxesOrigin->{-tplot/10,.0295},AxesLabel->{"Years",None},PlotStyle->{{Thickness[.007],Dashing[{.05,.0}]},{Thickness[.007],Dashing[{.02,0.02}]},{Thickness[.007],Dashing[{.01,0.01}]}},PlotLabel->"Foreign Currency Balances",DisplayFunction->Identity];
pk=Plot[{If[t<0,(Subscript[k, o]),k[t]]/.gminsol,(Subscript[k, f])},{t,-tplot/10,tplot},PlotRange->{.99,1.07},AxesOrigin->{-tplot/10,.99},AxesLabel->{"Years",None},PlotStyle->{{Thickness[.007],Dashing[{.05,.0}]},{Thickness[.007],Dashing[{.02,0.02}]},{Thickness[.007],Dashing[{.01,0.01}]}},PlotLabel->"Capital Stock",DisplayFunction->Identity];
pY=Plot[{If[t<0,Y[(Subscript[k, o])],Y[k[t]]]/.gminsol,Y[(Subscript[k, f])]},{t,-tplot/10,tplot},PlotRange->{.373,.388},AxesOrigin->{-tplot/10,.373},AxesLabel->{"Years",None},PlotStyle->{{Thickness[.007],Dashing[{.05,.0}]},{Thickness[.007],Dashing[{.02,0.02}]},{Thickness[.007],Dashing[{.01,0.01}]}},PlotLabel->"Real Output",DisplayFunction->Identity];
pr=Plot[{If[t<0,rfunc[(Subscript[m, o]),(Subscript[c, o]),(Subscript[F, o])],rfunc[m[t],c[t],F[t]]]/.gminsol,rfunc[(Subscript[m, f]),(Subscript[c, f]),(Subscript[F, f])]},{t,-tplot/10,tplot},PlotRange->{.098,.118},AxesOrigin->{-tplot/10,.098},AxesLabel->{"Years",None},PlotStyle->{{Thickness[.007],Dashing[{.05,.0}]},{Thickness[.007],Dashing[{.02,0.02}]},{Thickness[.007],Dashing[{.01,0.01}]}},PlotLabel->"Real Interest Rate",DisplayFunction->Identity];
ppi=Plot[{If[t<0,pifunc[(Subscript[m, o]),(Subscript[c, o]),(Subscript[F, o])],pifunc[m[t],c[t],F[t]]]/.gminsol,pifunc[(Subscript[m, f]),(Subscript[c, f]),(Subscript[F, f])]},{t,-tplot/10,tplot},PlotRange->{.20,.70},AxesOrigin->{-tplot/10,.20},AxesLabel->{"Years",None},PlotStyle->{{Thickness[.007],Dashing[{.05,0}]},{Thickness[.007],Dashing[{.02,0.02}]},{Thickness[.007],Dashing[{.01,0.01}]}},PlotLabel->"Inflation",DisplayFunction->Identity];
Show[GraphicsArray[{{pc,pin},{pm,pF},{pY,pk},{pr,ppi}}]]

(* Tabulate the final solutions. *)
Labels={"m","c","i","k","F","Y","z","v","r","pi"};InitialSteadyState={(Subscript[m, o]),(Subscript[c, o]),δ (Subscript[k, o]),(Subscript[k, o]),(Subscript[F, o]),(Subscript[Y, o]),(Subscript[z, o]),(Subscript[v, o]),ρ,(Subscript[π, o])}//N;FinalSteadyState={(Subscript[m, f]),(Subscript[c, f]),δ (Subscript[k, f]),(Subscript[k, f]),(Subscript[F, f]),(Subscript[Y, f]),(Subscript[z, f]),(Subscript[v, f]),ρ,(Subscript[π, f])}//N;
(Dat=Table[{tt,m[t],c[t],in[t],k[t],F[t],Y[k[t]],z[t],v[t],rfunc[m[t],c[t],F[t]],pifunc[m[t],c[t],F[t]]}/.tt->t/.gminsol,{t,0,T+0,1/5}])//MatrixForm;
(Dat=Join[{Join[{"Initial Steady State"},InitialSteadyState],Join[{"Final Steady State"},FinalSteadyState],Join[{Time},Labels]},Flatten[Dat,1]])//N//MatrixForm

The solve forward shoot 3D Package

(* Interface *)
BeginPackage["solve`"];

EndPackage[];

Off[General::spell];Off[General::spell1];Off[Unset::norep];Off[Unset::wrsym];Off[Precision::precsm];Off[InterpolatingFunction::dmval];

(* Usage messages for the exported functions and the context itself *)
solve::usage = "solve[sys, {t}, vars, varso, varsf, varsendo, sysexpand, varspecial, init, pini] solves a system of first-order differential equations. \n
sys   = {xdot, ydot, ...} lists the right hand side of the system of first-order differential equations (with jump variables in the core dynamic system in beginning and state variables in the end).\n
vars  = {x, y, ...} is the list of variables in the system in the same order as the differential equations in sys.\n
varso = {(Subscript[x, o]), (Subscript[y, o]), ...} is the list of original steady values of variables in the system in the same order as vars.\n
varsf = {(Subscript[x, f]), (Subscript[y, f]), ...} is the list of final steady values of variables in the system in the same order as vars.\n
varsendo = {{w, wfunc, {x, y, ...}}, ...} is the list of functions {wfunc, ...} for the endogenous variables {w, ...} where wfunc is function of {x, y, ...} and {x, y, ...} are in the order in which they appear as arguments of wfunc. If no endogenous variables in the system, enter {{,{}}}.\n
sysexpand  = {funcj, {x, y, ...}} contains the name of function funcj that takes values of variables in core dynamic 
	system {x, y, ...} at time 0 as argument and returns a list containing values of other jump variables (in the same order as they appear in vars) that have been added to core dynamic system to obtain the expanded system sys. If core dynamic system has not been expanded, enter {}.\n
varspecial  = {v, ...} is the list of variables that can potentially take non-positive values. It is assumed that other variables in vars have been defined so that they only take positive values. If no variable can potentially take non-positive values, enter {}.\n
init  = {T, tol, precision, epsilon} sets the required duration of the solution to T, the required accuracy of solution to tol, starting precision to precision, and radius of initial ϵ-ball to epsilon. Note that setting T larger than necessary can significantly increase computation time. \n
pini  = {a, b, c} gives center of the initial ring in the jump space. The typical value is {(Subscript[x, o]), (Subscript[y, o]), (Subscript[z, o])}, the initial steady state values of the jump variables in the core dynamic system. A better guess may somewhat help reduce computational time but, typically, is not necessary.";

(* Error messages for the exported objects *)
{  solve::sysDESys = "You called `1` with first argument `2` (the system of differential equations) that has variables other than those defined in `3` and `4`.";

solve::sysLength = "You called `1` with lengths of `2`, `3`, `4` and `5` not equal.";

solve::vars = "You called `1` with third argument `2` that is not a list of variables.";

solve::origVals = "You called `1` with fourth argument `2` that is not a list of numbers.";

solve::finVals = "You called `1` with fifth argument `2` that is not a list of numbers.";

solve::prec = "You called `1` with fourth and fifth arguments `2` and `3` that have numbers with precision less than 100.";

solve::endoVars = "You called `1` with sixth argument `2` that is not a list of lists of the form {w, wfunc, {x, y, ...}} where w is an endogenous variable and {x, y, ...} are the variables in the system as defined in `3`.";

solve::initVals = "You called `1` with last argument `2` that is not a list of numbers.";

solve::sysError = "You called `1` with a system that does not have three negative eigenvalues when linearized."};

(* Implementation *)
Begin["solve`Private`"];

(* Definition of the exported functions *)
solve[asys_List,{at_},avars_List,avarso_List,avarsf_List,avarsendo_List,asysexpand_List,avarspecial_List,ainit_List,apini_List]:=Block[{sys,t,vars,varso,varsf,varsendo,init,sysexpand,varspecial,T,tol,precision,delta,Tmax,firstsol,precisionini,nringsini,system,boxvars,box,A,localmin,pmin,epsilon,digits,ntrajini,gmindis,gtmin,gminsol,initsol,ini,dis,tsol,tdis,doneflag},

{sys,{t},vars,varso,varsf,varsendo,sysexpand,varspecial,init,pmin}={asys,{at},avars,avarso,avarsf,avarsendo,asysexpand,avarspecial,ainit,apini};
{T,tol,precision,epsilon}=init;Tmax=4000;firstsol = 1;precisionini=precision;nringsini=12;

(* Check input *)
checkInput[];

(* Make the system from input data and create box over which to solve the system *)
system=makeSys[];
boxvars=Map[Part[#,Flatten[Position[Map[MemberQ[varspecial,#]&,vars],False]]]&,{vars,varso,varsf}];
box:=Apply[Or,Join[MapThread[First[Map[(#1),{t}]]<((Min[#2,#3])/20)&,boxvars],MapThread[First[Map[(#1),{t}]]>5 Max[#2,#3]&,boxvars]]];

(* SOLVE FOR THE GLOBAL NONLINEAR SADDLEPATH *);
While[True,
localmin=0;

Label[begin];
$MinPrecision=digits=precision;ntrajini=11;gmindis=tol;
If[firstsol==0,{pmin,epsilon}=(<<solution2)[[{4,5}]]];

(* Calculate the minimum distance for the trajectory of the nonlinear system that starts from the center of the initial ϵ-ring. *)
{gmindis,gtmin,gminsol,pmin}=compdis[pmin,0,0,0];

Print["AT THE CENTER OF THE INITIAL BALL: ","    Minimum Distance from Final Steady State = ",gmindis//N,"    Time of Closest Approach = ",gtmin//N];

(* SOLVE - Solve for the global nonlinear saddle path *);
If[localmin==1,localmin=0;pmin=pmin+epsilon/10{1/2,1/2}];

{gminsol,gtmin,gmindis,pmin,epsilon}=compdis3[pmin,epsilon];

If[gtmin==(10^100),
Print["Either the system is highly ill-conditioned or you might have encountered a local minima. Perturing the center of the initial ring and resetting 'precision' to its initial value"];precision=precisionini;localmin=1;Goto[begin]];

If[doneflag==1,Goto[end]];

];
 
(* DURATION CHECK - Check the duration for which the solution is accurate within tolerance 'tol'.*);
Label[end];
durationCheck[];

{gtmin,gminsol}];

checkInput[___]:=Block[{},

If[Not[Length[sys]==Length[vars]==Length[varso]==Length[varsf]],
	Message[solve::sysLength,"'solve'","'sys'","'vars'","'varso'","'varsf'"];Abort[]];

If[Not[Length[vars]==Length[Cases[vars,_Symbol]]],Message[solve::vars,"'solve'","'vars'"];Abort[]];

If[Not[Length[varso]==Length[Cases[varso,_?NumericQ]]],Message[solve::origVals,"'solve'","'varso'"];Abort[]];

If[Not[Length[varsf]==Length[Cases[varsf,_?NumericQ]]],Message[solve::finVals,"'solve'","'varsf'"];Abort[]];

If[Min[Map[Precision[#]&,Join[varso,varsf]]]<100,Message[solve::prec,"'solve'","'varso'","'varsf'"];Abort[]];

If[Not[Length[varsendo]==Length[Cases[varsendo,{_Symbol,_Symbol,_List}]]]||Not[Length[Flatten[varsendo]]==Length[Cases[Flatten[varsendo],_Symbol]]],
	Message[solve::endoVars,"'solve'","'varsendo'","'vars'"];Abort[];];

If[Not[Length[init]==Length[Cases[init,_?NumericQ]]],Message[solve::initVals,"'solve'","'init'"];Abort[]];

];

makeSys[___]:=Block[{system},

system=sys/. Thread[MapThread[(#1)&,Transpose[varsendo]]->
					MapThread[ToExpression[StringJoin[ToString[ #2],"[",ToString[#3],"]"]] /.List->Sequence&, Transpose[varsendo]],
						List];

If[Length[Cases[system/.t->0/.Thread[vars->varsf,List],_?NumberQ]]≠Length[vars],
	Message[solve::sysDESys,"'solve'","'sys'","'vars'","'varsendo'"];Abort[];];

system=Thread[Map[ToExpression[StringJoin[ToString[#],"'","[",ToString[t],"]"]]&,vars]==  
			(system=system/.Thread[vars->Map[ToExpression[StringJoin[ToString[#],"[",ToString[t],"]"]]&,vars],List]),
				List];
system];

distance[sol_,tsol_,tini_]:=Block[{dis,dislast,tdis,steps,tbeg,tend,intsz,traj,min,norm},
dis= (10^7);dislast=(10^8);If[firstit==1,steps=Max[4Floor[tsol+1],450] ,steps=Max[4Floor[tsol+1-tini],50] ];
tbeg=tini;tend=tsol;intsz=tend-tbeg;
norm=Table[vars[[i]][t]-varsf[[i]],{i,1,Length[vars]}];

While[ (dislast-dis)>gmindis (10^(-7)),
	dislast=dis;
	traj=Table[Norm[norm/.t->s/.sol,2],{s,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}];

(* Set up for solving for the global nonlinear saddle path *)
compdis3[p_,ϵ_]:=Block[{bmin,btmin,bminsol,bpmin,bit,shrink,difftmin,notsatistmin,satistmin,diffprec,notsatisprec,satisprec,increaseprec,congratulations,done,step2,step3,iflag,rminpos,ringtrack},
iflag=0;bit=1;shrink=False;rminpos={};ringtrack=0;

(* Messages *)
difftmin:=Print["Difference in values of variables between two paths with (Subscript[t, min])s differing by at least 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! The variables changed by less than 'tol' due to latest increase in 'precision'.";
done = "You are done! Hope to see you in another model. Continue with the next cell.";

If[firstsol==0,bit=0;Goto[errorcheck]];

Label[begin];
Print[""];Print["BEGINNING SEARCH ON BALL ",bit," WITH CENTER ",N[pmin,precision]];
precision=precisionini;$MinPrecision=digits=precision;

{bmin,btmin,bminsol,bpmin}=compdis2[pmin,epsilon,0,-(π/2) ];

If[btmin==(10^100)&& precision≤80,
	precision=precision+4;$MinPrecision=digits=precision;Print["Increased 'precision' by 4 digits to ",precision," digits."];Goto[begin],
	If[btmin==(10^100),gtmin=btmin;Goto[end]; ]];

If[bmin≤gmindis,
	gmindis=bmin;gtmin=btmin;gminsol=bminsol;pmin=bpmin;
	If[shrink==True,epsilon=epsilon/2,epsilon=epsilon/1],
	If[bit≠1,shrink=True];
	epsilon=epsilon/2];

Print["AT THE END OF SEARCH ON BALL ",bit,":    CURRENT MINIMUM DISTANCE = ",gmindis//N,"     TIME OF CLOSEST APPROACH = ",gtmin];

(* ERROR CHECK - Check for the accuracy of the solution over the time period of interest. *) 
Label[errorcheck];
If[firstsol==1,
	If[gtmin>T+5 && iflag==0,iflag=1;{gminsol,gtmin,gmindis,pmin,5epsilon}>>solution1;Print[""];Print[step2];];
	If[gtmin>Max[T+5 ,(<<solution1)[[2]]+5]&& iflag==1 ,
		{gminsol1,gtmin1}=(<<solution1)[[{1,2}]];
		If[Max[sysDiff[0,T]]>tol,
			Print[""];difftmin;Print[notsatistmin];
			{gminsol,gtmin,gmindis,pmin,5epsilon}>>solution1,
			
			Print[""];difftmin;Print[satistmin];
			maxErrorUptoT=Max[sysDiff[0,T]];tempT=IntegerPart[gtmin1];
			For[j=T+1,j≤IntegerPart[gtmin1],j++,If[Max[sysDiff[j-1,j]]>tol,tempT=j-1;Break[]];];
			Print[""];Print[step3];firstsol=0;precision=precision+4;increaseprec;
			{gminsol,gtmin,gmindis,pmin,5epsilon}>>solution2;doneflag=0(*THIS DONEFLAG COMMAND SKIPS STEP 3 if set to 1*);
			Goto[end];]],
	
	If[gtmin>((<<solution2)[[2]]-1)&& gmindis<(101/100)(<<solution2)[[3]],
		gminsol1=(<<solution2)[[1]];
		If[Max[sysDiff[0,T]]>tol,
			Print[""];diffprec;Print[notsatisprec];precision=precision+4;increaseprec,
			Print[""];diffprec;Print[satisprec];Print[""];Print[congratulations];Print[done];doneflag=1];
			Goto[end];];
		(* solution1 stores the final solution from step 3 and solution2 stores the final solution from step 2 unless there is a failure of precision in step 3. *);
];

bit++;Goto[begin];

Label[end];
{gminsol,gtmin,gmindis,pmin,5epsilon}];

Clear[compdis2];
compdis2[p_,epsilon_,thetabeg_,phibeg_]:=Block[{rmin,rmin1,rmin2,rmindis,nrings,jr,ir,jr1,jr2,sol,dis,tdis,rtmin,jrmin,rminsol,rtmin2,rtab,jrmin2,rminsol2,rmindis2,theta,jrdis,jrtdis,jrdissol,jrp,rp,trajtrack,minpos,firstringit,},

nrings=nringsini;firstringit=1;minpos={};trajtrack=0;rp=p;

solvering:=(Label[restart];
	{dis,tdis,sol,rmin}=compdis[rp,epsilon,theta,phibeg];
	If[tdis==(10^100)&& precision≤80,
	precision=precision+4;$MinPrecision=digits=precision;If[firstringit==1,precisionini=precision];
	Print["Increased 'precision' by 4 digits to ",precision," digits."];rp=SetPrecision[rp,precision];Goto[restart],
	If[tdis==(10^100),rtmin2=tdis;Goto[end1] ]]);

For[jr=1,jr≤ nrings ,jr++,
	If[ringtrack==0||(ringtrack==1 && Apply[Or,Map[jr==# &,rminpos]]),
	theta= thetabeg+((2π jr)/nrings);
	solvering;
	jrdis[jr]=dis;jrtdis[jr]=tdis;jrdissol[jr]=sol;jrp[jr]=rmin,jrdis[jr]=(10^7)];
];

firstringit=0;
If[Length[rminpos]==2 && iflag==1,
ringtrack= 1,
rtab=Table[jrdis[Mod[ir-1,nrings]+1]-jrdis[Mod[ir,nrings]+1],{ir,1,nrings}];rminpos=Flatten[Position[Table[rtab[[Mod[ir-2,nrings]+1]]>0 &&rtab[[Mod[ir-1,nrings]+1]]<0,{ir,1,nrings}],True]];Print[" "];];

{rmindis,rmindis2}={(10^7),(10^7)};

For[ir=1,ir≤Length[rminpos],ir++,
	Print["Initial Ring = ",ir,"   {rmindis, rtdis, rjmin, nrings} = ",{jrdis[rminpos[[ir]]],jrtdis[rminpos[[ir]]]//N,rminpos[[ir]],nrings}]];Print[" "];

For[ir=1,ir≤Length[rminpos],ir++,
{rmindis,rtmin,jrmin,nrings,rminsol}={jrdis[rminpos[[ir]]],jrtdis[rminpos[[ir]]],rminpos[[ir]],nringsini,jrdissol[rminpos[[ir]]]};
rmin1=jrp[rminpos[[ir]]];
minpos={};trajtrack=0;precision=precisionini;$MinPrecision=digits=precision;

While[True,
{jr1,jr2}=2{jrmin-1,jrmin+1};{nrings,jrmin}=2{nrings,jrmin};
For[jr=jr1+1,jr≤ jr2  ,jr+=2,
	theta= thetabeg+((2π jr)/nrings);
	solvering;
	If[ Abs[rmindis-dis]≤(gmindis/100000), Goto[end]; ];
	If[dis≤rmindis,
		rmindis=dis;rtmin=tdis;jrmin=jr;rminsol=sol;rmin1=rmin;If[bit≥5&&rmindis<(75/100)gmindis,Goto[end];];Break[];];
];
];

Label[end];
If[rmindis<rmindis2,rmindis2=rmindis;rtmin2=rtmin;jrmin2=jrmin;rminsol2=rminsol;rmin2=rmin1];

If[bit ≥5 && rmindis2<gmindis,Break[]];
];Print[""];

Label[end1];
Print["Ball = ",bit,"    Minimum Distance = ",rmindis2//N,"    Minimum point = ",N[rmin2,precision],"    No. of Subdivisions of ϵ-ball > ",jrmin2];{rmindis2,rtmin2,rminsol2,rmin2}];

(* Set up for solving for the global nonlinear saddle path *)
solvesys[ini_,t1_,t2_]:=Block[{sol},
	sol=NDSolve[{system,ini},vars,{t,t1,t2},Method->{"EventLocator", "Event"->box,"EventAction":> Throw[Null, "StopIntegration"],"EventLocationMethod"->"StepBegin"},
		MaxSteps->Infinity,WorkingPrecision->digits,PrecisionGoal->Min[10,digits/2],AccuracyGoal->Infinity];

sol];
compdis[p_,ϵ_,theta_,phibeg_]:=Block[{solveit,phi,ntraj,pini,pexp,ini,sol,dis,tdis,j,jdis,jtdis,dissol,jp,tab,mindis,mindis2,j1,j2,ii,tmin,jmin,minsol,tmin2,pmin1,jmin2,minsol2,pmin2},

ntraj=ntrajini;firstit=1;

solveit:=(phi = phibeg+((π )/ntraj)j;pini=p+ϵ{1Cos[theta]Cos[phi],1Sin[theta]Cos[phi],1Sin[phi]};
If[Length[sysexpand]==0,pexp={},pexp=First[sysexpand][Apply[Sequence,Last[sysexpand]]];];

(* IF CLAUSE FOR OUT OF BOX SO THAT IT IS DONE ONLY FOR SOME INITIAL BALLS *)If[bit≤5,Table[vars[[i]][0]=.,{i,1,Length[vars]}];Evaluate[Table[vars[[i]][0],{i,1,Length[vars]}]]=Apply[Sequence,{pini,pexp,Drop[varso,3+Length[pexp]]},1];
	If[(box/.t->0)== True,Print["Out of Box."];{dis,tdis}={(10^7),100};Goto[skip]];Table[vars[[i]][0]=.,{i,1,Length[vars]}];];

ini=Thread[Table[vars[[i]][0],{i,1,Length[vars]}]==Apply[Sequence,{pini,pexp,Drop[varso,3+Length[pexp]]},1],List];
sol=solvesys[ini,0,Tmax];{dis,tdis}=distance[sol,sol[[1,1,2,1,1,2]],If[firstit==1,0,Max[0,Min[tmin/2,tmin-5]]]];Label[skip]);

If[ϵ==0,solveit;mindis2=dis;tmin2=tdis;pmin2=p;minsol2=sol;Goto[end2]];

For[j=-1,j≤ ntraj+1  ,j++,
	If[trajtrack==0||(trajtrack==1 && Apply[Or,Map[j==# &,minpos-1]]),
		solveit;jdis[j]=dis;jtdis[j]=tdis;dissol[j]=sol;jp[j]=pini,jdis[j]=(10^7)];
];

firstit=0;
If[jr≥ 5nringsini,
trajtrack= 1,

tab=Table[jdis[ii-1]-jdis[ii],{ii,0,ntraj+1}];
minpos=Flatten[Position[Table[tab[[ii]]>0 &&tab[[ii+1]]<0,{ii,1,ntraj}],True]];];

Print["Local minima from initial search over the ring."];
For[ii=1,
		ii≤Length[minpos],ii++,
		Print[ii,":" ,"   {mindis, tdis, jmin, ntraj} = ",{jdis[minpos[[ii]]-1],jtdis[minpos[[ii]]-1]//N,minpos[[ii]]-1,ntraj}];
];

mindis=mindis2=(10^7);j1=1;j2= ntraj;

For[ii=1,ii≤Length[minpos],ii++,
	{mindis,tmin,jmin,ntraj,minsol}={jdis[minpos[[ii]]-1],jtdis[minpos[[ii]]-1],minpos[[ii]]-1,ntrajini,dissol[minpos[[ii]]-1]};
	pmin1=jp[minpos[[ii]]-1];

	While[True,
	{j1,j2}=2{jmin-1,jmin+1};{ntraj,jmin}=2{ntraj,jmin};
	For[j=j1+1,j≤ j2  ,j+=2,
		solveit;
		If[ Max[Abs[Apply[(#1/#2-1)&,{pini,pmin1}]]]<(10^(-precision)),
			Print["You may be losing precision in calculations."] ;
			tmin2=(10^100);Goto[end2];];
		If[ Abs[mindis-dis]≤gmindis/100000,Goto[end1];];
		If[dis≤mindis,mindis=dis;tmin=tdis;jmin=j;minsol=sol;pmin1=pini;Break[];];
	];
	];

	Label[end1];

	If[mindis<mindis2,
	mindis2=mindis;tmin2=tmin;jmin2=jmin;minsol2=minsol;pmin2=pmin1];
	If[mindis2<gmindis ,Break[]];

];
Print["Ring ",jr,":    Minimum Distance = ",mindis2//N,"    Time of Closest Approach = ",tmin2//N,"    No. of Subdivisions of ϵ-ring > ",jmin2];

Label[end2];
{mindis2,tmin2,minsol2,pmin2}];

sysDiff[t1_?NumberQ,t2_?NumberQ]:=Table[Max[Abs[Table[((vars[[i]][t]/.gminsol)-(vars[[i]][t]/.gminsol1))[[1]],{t,t1,t2,1/10}]]],{i,1,Length[vars]}];

durationCheck[___]:=Block[{j},

gminsol1=(<<solution2)[[{1}]]; {gminsol,gtmin}=(<<solution1)[[{1,2}]];

For[j=T+1,j≤IntegerPart[tempT],j++,
	If[Max[maxErrorUptoT,Max[sysDiff[j-1,j]]]>tol,
		Print["The solution is accurate within tolerance 'tol' = ",tol//N," upto ",j-1," years."];Break[]];];

If[j==IntegerPart[tempT]+1,Print["The solution is accurate within tolerance 'tol' = ",tol//N," upto ",j-1," years. Continue with next cell."];];
Print["The solution upto ",T," years is accurate to within ",Max[maxErrorUptoT,Max[sysDiff[0,T]]]//N,". Continue with next cell."];

];

End[];

(* Epilog *)
Protect[Evaluate[$Context<>"*"]];