This file implements the algorithm 'Forward Shoot 2D'  outlined in paper titled, "Smart Forward Shooting" which solves models with two 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 2 of the paper.

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

(* 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),σ=(1/3),ρ=(1/10),β=(1/2),δ=(1/20)}//N                                   (* Parameters of the model *)

{(Subscript[bbyY, o])=(25/100),(Subscript[mbyY, 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 *)

Steady State

(* Define utility function and production function *)
U[c_,m_]:=((1 )/(1-(1/τ)))(((((((c^(1-(1/σ)))+g (m^(1-(1/σ)))))^(σ/(σ-1)))))^(1-(1/τ)));(Subscript[U, 1])[c_,m_]:=Evaluate[D[U[c,m],c]];(Subscript[U, 2])[c_,m_]:=Evaluate[D[U[c,m],m]];
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]),b=(Subscript[bbyY, o]) (Subscript[Y, o]),(Subscript[z, o])=(Subscript[π, o])(Subscript[m, o])+(Subscript[v, o])Y[(Subscript[k, o])]-ρ b};g=((Subscript[π, o])+ρ)(((((Subscript[m, o])/(Subscript[c, o]))))^(1/σ));{(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[k, o]))/(Subscript[Y, o]))},{(Subscript[a, 0]),1},{(Subscript[a, 1]),1},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])}=.;{(Subscript[m, f]),(Subscript[π, f])}={(Subscript[m, f]),(Subscript[π, f])}/.FindRoot[{g(((((Subscript[m, f])/(Subscript[c, f]))))^(-1/σ))==(Subscript[π, 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]),(Subscript[π, o])},AccuracyGoal->digits,WorkingPrecision->digits,MaxIterations->1000];
{{"m","c","k","Y","b","z","pi","c/m","m/Y"},{(Subscript[m, 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[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

(* 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; 
kdot = Y[k]-c-δ k;
λdot =(ρ-((1-v)(Subscript[Y, 1])[k]-δ))λ;
cdot=dcin/.Solve[Dt[(Subscript[U, 1])[c,m]==λ]/.{Dt[c]->dcin,Dt[m]->mdot,Dt[k]->kdot,Dt[v]->vdot,Dt[λ]->λdot},dcin][[1]]/.λ->(Subscript[U, 1])[c,m];
system  ={mdot,cdot,kdot};
};

(* 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[k_,t_]:=(1-v[t])(Subscript[Y, 1])[k]-δ;
pifunc[m_,c_, k_,t_]:=g((((m/c)))^(-(1/σ)))-rfunc[k,t];
};

(* 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]);

(* INITIALIZATION CELL - Initialization for finding the required trajectory. *)
{  
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 ϵ-ring*)
pini={(Subscript[m, o]),(Subscript[c, o])};                                    (* center of the initial ϵ-ring *);
};
(* 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,k},{(Subscript[m, o]),(Subscript[c, o]),(Subscript[k, o])},{(Subscript[m, f]),(Subscript[c, f]),(Subscript[k, f])},{{r,rfunc,{k,t}},{pi,pifunc,{m,c,k,t}},{z,z,{t}},{v,v,{t}}},{},{},{T,tol,precision,epsilon},pini]
Print["Maximum amount of memory used in computing the solution: ",MaxMemoryUsed[]/1000000MB//N]

(* Plot the final solutions. *)
tplot= T+5;
pc=Plot[{If[t<0,(Subscript[c, o]),c[t]]/.gminsol,(Subscript[c, f])},{t,-tplot/10,tplot},PlotRange->{.245,.262},AxesOrigin->{-tplot/10,.245},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];
pm=Plot[{If[t<0,(Subscript[m, o]),m[t]]/.gminsol,(Subscript[m, f])},{t,-tplot/10,tplot},PlotRange->{.0188,.026},AxesOrigin->{-tplot/10,.0188},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];
pk=Plot[{If[t<0,(Subscript[k, o]),k[t]]/.gminsol,(Subscript[k, f])},{t,-tplot/10,tplot},PlotRange->{.99,1.09},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->{.298,.316},AxesOrigin->{-tplot/10,.298},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,r[(Subscript[k, o]),0],rfunc[k[t],t]]/.gminsol,rfunc[(Subscript[k, f]),Infinity]},{t,-tplot/10,tplot},PlotRange->{.098,.106},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[k, o]),0],pifunc[m[t],c[t],k[t],t]]/.gminsol,pifunc[(Subscript[m, f]),(Subscript[c, f]),(Subscript[k, f]),Infinity]},{t,-tplot/10,tplot},PlotRange->{.20,.57},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];

If[$VersionNumber==5.2,
Show[GraphicsArray[{{pc,pm},{pY,pk},{pr,ppi}}]],
Show[GraphicsGrid[{{pc,pm},{pY,pk},{pr,ppi}}]]]

Labels={"m","c","i","k","Y","z","r","pi"};InitialSteadyState={(Subscript[m, o]),(Subscript[c, o]),δ (Subscript[k, o]),(Subscript[k, o]),(Subscript[Y, o]),(Subscript[z, o]),ρ,(Subscript[π, o])}//N;FinalSteadyState={(Subscript[m, f]),(Subscript[c, f]),δ (Subscript[k, f]),(Subscript[k, f]),(Subscript[Y, f]),(Subscript[z, f]),ρ,(Subscript[π, f])}//N;
(Dat=Table[{tt,m[t],c[t],Y[k[t]]-c[t],k[t],z[t],Y[k[t]],rfunc[k[t],t],pifunc[m[t],c[t],k[t],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 2D 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 ϵ-ring to epsilon. Note that setting T larger than necessary can significantly increase computation time. \n
pini  = {a, b} gives center of the initial ring in the jump space. The typical value is {(Subscript[x, o]), (Subscript[y, 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,Tmax,firstsol,precisionini,system,boxvars,box,localmin,pmin,epsilon,digits,ntrajini,gmindis,gtmin,gminsol,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;

(* Check the input from the user for errors. *)
checkInput[];

(* Make the system from the input data. *)
system=makeSys[];

(* Create box over which to solve the system. *) 
(* If during computation, a path goes out of this box, its further computation is terminated. *) 
(* The box is formed with all the variables except those that can potentially take non-positive values. *) 
boxvars=Map[Part[#,Flatten[Position[Map[MemberQ[varspecial,#]&,vars],False]]]&,{vars,varso,varsf}];
box:=Evaluate[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=20;gmindis=tol;
If[firstsol==0,{pmin,epsilon}=(<<solution2)[[{4,5}]]];

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

Print["AT THE CENTER OF THE INITIAL RING: ","    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}=compdis2[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[{},

(* Check if number of differential equations in the system, the number of variables in the system, and the lists containing their initial and final steady state values have consistent dimensions. *)
If[Not[Length[sys]==Length[vars]==Length[varso]==Length[varsf]],
	Message[solve::sysLength,"'solve'","'sys'","'vars'","'varso'","'varsf'"];Abort[]];

(* Check if the variable names defined are valid Mathematica symbols. *)
If[Not[Length[vars]==Length[Cases[vars,_Symbol]]],Message[solve::vars,"'solve'","'vars'"];Abort[]];

(* Check if the initial steady state values are numeric. *)
If[Not[Length[varso]==Length[Cases[varso,_?NumericQ]]],Message[solve::origVals,"'solve'","'varso'"];Abort[]];

(* Check if the initial steady state values are numeric. *)
If[Not[Length[varsf]==Length[Cases[varsf,_?NumericQ]]],Message[solve::finVals,"'solve'","'varsf'"];Abort[]];

(* Check if the initial and final steady state values have precision of 100. *)
If[Min[Map[Precision[#]&,Join[varso,varsf]]]<100,Message[solve::prec,"'solve'","'varso'","'varsf'"];Abort[]];

(* Check if the information about endogenous variables in the system is provided in the required format. *)
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[];];

(* Check if the initialization information has only numeric quantities. *)
If[Not[Length[init]==Length[Cases[init,_?NumericQ]]],Message[solve::initVals,"'solve'","'init'"];Abort[]];

];

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

(* Replace the endogenous variables in the system with the corresponding user defined functions. *)
system=sys/. Thread[MapThread[(#1)&,Transpose[varsendo]]->
					MapThread[ToExpression[StringJoin[ToString[ #2],"[",ToString[#3],"]"]] /.List->Sequence&, Transpose[varsendo]],
						List];

(* Check if the system contains variables other than those specified by the user. *)
If[Length[Cases[system/.t->0/.Thread[vars->varsf,List],_?NumberQ]]≠Length[vars],
	Message[solve::sysDESys,"'solve'","'sys'","'vars'","'varsendo'"];Abort[]s;];

(* Create the final system of equation with time subscript 't' added to all the variables in the system. *)
system=Thread[Map[ToExpression[StringJoin[ToString[#],"'","[",ToString[t],"]"]]&,vars]==  
			(system=system/.Thread[vars->Map[ToExpression[StringJoin[ToString[#],"[",ToString[t],"]"]]&,vars],List]),
				List];

system];

(* This function computes the minimum distance of a trajectory from the final steady state over interval [tini, tsol]. *)
distance[sol_,tsol_,tini_]:=Block[{dis,dislast,tdis,steps,tbeg,tend,intsz,traj,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;

(* Define the Euclidean norm for computing distance. *)
norm=Table[vars[[i]][t]-varsf[[i]],{i,1,Length[vars]}];

(* Use repeated subdivision of search interval to compute the minimum. *)
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 *)
compdis2[p_,ϵ_]:=Block[{rmin,rtmin,rminsol,rpmin,difftmin,notsatistmin,satistmin,diffprec,notsatisprec,satisprec,increaseprec,congratulations,done,step2,step3,rit,iflag,shrink,pmin,epsilon},

iflag=0;trajtrack=0;rit=1;shrink=False;{pmin,epsilon}={p,ϵ};

(* 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,rit=0;Goto[errorcheck]];

Label[begin];
Print[""];Print["BEGINNING SEARCH ON RING ",rit," WITH CENTER ",N[pmin,precision]];

(* Find the minimum of the distance mapping on the current ring. *)
{rmin,rtmin,rminsol,rpmin}=compdis[pmin,epsilon];

(* Check if there was a failure of precision in finding the minimum. If so increase the precision and try again. *) (* If precision goes above 80 exit compdis2 as either system is too ill-conditioned or you encountered a local minima. *)
If[rtmin==(10^100)&& precision≤80,
	precision=precision+4;$MinPrecision=digits=precision;Print["Increased 'precision' by 4 digits to ",precision," digits."];Goto[begin],
	If[rtmin==(10^100),gtmin=rtmin;Goto[end] ]];

(* Compare the minimum for the current ring with the previous minimum. *)
(* If no improvement, then shrink the ring but do not change the center of the ring for the next search. *)(* If current minimum is smaller, make it the center for the next ring. Also if shrink is True, shrink the ring. *)
If[rmin≤gmindis,
	gmindis=rmin;gtmin=rtmin;gminsol=rminsol;pmin=rpmin;
	If[shrink==True,epsilon=epsilon/2,epsilon=epsilon/1],
	If[rit≠1,shrink=True];
	epsilon=epsilon/2];

Print["AT THE END OF SEARCH ON RING ",rit,":    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, a candidate solution has been obtained. Proceed to step 2. *)
	If[gtmin>T+5 && iflag==0,iflag=1;{gminsol,gtmin,gmindis,pmin,5epsilon}>>solution1;Print[""];Print[step2];];
	
	(* If found another solution with gtmin exceeding previous gtmin by 5 years,... *)
       (* ... a solution has been obtained for checking the accuracy of the candidate solution. *) 
	If[gtmin>Max[T+5 ,(<<solution1)[[2]]+5]&& iflag==1 ,
		{gminsol1,gtmin1}=(<<solution1)[[{1,2}]];
		(* Check if accuracy check is passed.*)
		If[Max[sysDiff[0,T]]>tol,
			(* If not then, the current solution becomes the new candidate solution. *)
			(* In this case, look for another solution whose duration exceeds the current candidate solution by 5 years.*)
			Print[""];difftmin;Print[notsatistmin];
			{gminsol,gtmin,gmindis,pmin,5epsilon}>>solution1,
			
			(* If accuracy check is passed, proceed to step 3. *)
			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;
			Goto[end]]],

	(* Check if found a solution with similar or better duration and minimum distance with increased precision in step 3. *)
	If[gtmin>((<<solution2)[[2]]-1)&& gmindis<(101/100)(<<solution2)[[3]],
		gminsol1=(<<solution2)[[1]];
		(* Check if precision check is passed.*)
		If[Max[sysDiff[0,T]]>tol,
			(* If not then, the current solution becomes the new candidate solution. *)
			{gminsol,gtmin,gmindis,pmin,epsilon}>>solution2;
			Print[""];diffprec;Print[notsatisprec];precision=precision+4;increaseprec,

			(* If precision check is passed, you are done. *)
			{gminsol,gtmin,gmindis,pmin,epsilon}>>solution1;
			Print[""];diffprec;Print[satisprec];Print[""];Print[congratulations];Print[done];doneflag=1];
			Goto[end]];

];

rit++;Goto[begin];

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

(* Set up for solving for the global nonlinear saddle path. *)

(* The function solvesys solves the system for period [t1,t2] with initial conditions ini at t = t1. *)
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];

(* The function compdis identifies the point or initial conditions on the current ring in the jump space which yields the trajectory that comes closest to the final steady state. This constitutes the basic unit of circle search. *)
compdis[p_,ϵ_]:=Block[{solveit,theta,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 is called at each point on the current ring to solve for the corresponding trajectory. *)
solveit:=(
(* compute the exact location of the current point on the ring. *)
theta = 2π j/ntraj;pini=p+ϵ{Sin[theta], Cos[theta]};

(* If the core dynamic system has been expanded, generate the initial values of the variables added. *)
If[Length[sysexpand]==0,pexp={},pexp=First[sysexpand][Apply[Sequence,Last[sysexpand]]];];

(* Check if the starting point of a trajectory is itself out of the box. In that case, do not solve for the trajectory. This check is done only for some initial rings. *)If[rit≤5,Table[vars[[i]][0]=.,{i,1,Length[vars]}];Evaluate[Table[vars[[i]][0],{i,1,Length[vars]}]]=Apply[Sequence,{pini,pexp,Drop[varso,2+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]}];];

(* Set up the initial conditions corresponding to the current point on the ring. *)
ini=Thread[Table[vars[[i]][0],{i,1,Length[vars]}]==Apply[Sequence,{pini,pexp,Drop[varso,2+Length[pexp]]},1],List];

(* Call solvesys to solve for the required trajectory and calculate its minimum distance from the final steady state. *)
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]);

(* Execute this command, and only this command, if compdis is called from solve for a zero radius ring. *)
If[ϵ==0,solveit;mindis2=dis;tmin2=tdis;pmin2=p;minsol2=sol;Goto[end2]];

(* First pass on the ring to locate the minimum of the distance mapping on the ring. *)
For[j=1,j≤ ntraj  ,j++,
	If[trajtrack==0||(trajtrack==1 && (j==minpos[[1]]||j==minpos[[2]])),
		solveit;jdis[j]=dis;jtdis[j]=tdis;dissol[j]=sol;jp[j]=pini,jdis[j]=(10^7)];
];

firstit=0;

(* Determine the location of minima on the ring from the initial pass on the ring *)
(* If there are exactly two, then in future computations, initial search will be limited to those trajectories only. *)
If[iflag==1 &&Length[minpos]==2,
trajtrack= 1,

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

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

(* Now implement bisection algorithm around each local minima identified to find the minimum on the ring. *)
For[ii=1,ii≤Length[minpos],ii++,
	(* Get the information about current local minima. *)
	{mindis,tmin,jmin,ntraj,minsol}={jdis[minpos[[ii]]],jtdis[minpos[[ii]]],minpos[[ii]],ntrajini,dissol[minpos[[ii]]]};
	pmin1=jp[minpos[[ii]]];
	
	(* Do repeated bisection until minimum of the distance mapping does not change by more that (10^(-5))gmindis. *)
	While[True,
	{j1,j2}=2{jmin-1,jmin+1};{ntraj,jmin}=2{ntraj,jmin};
	
	For[j=j1+1,j≤ j2  ,j+=2,
		(* Call solveit to solve for a trajectory. *)
		solveit;
		(* If points on the ring become too close (distance < (10^(-preicision)))there is risk of loss of precision. Exit compdis. *)
		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 ",rit,":    Minimum Distance = ",mindis2//N,"    Time of Closest Approach = ",tmin2//N,"    No. of Subdivisions of ϵ-ring > ",jmin2];

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

(* This function checks the maximum absolute difference between different variables in the system for two candidate solutions. *)
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]}];

(* This function gives the duration for which the finally accepted solution is accurate within specified tolerance. It also gives the actual tolerance achieved for the specified time duration, T. *)
durationCheck[___]:=Block[{j},
(* Get the two paths that were compared during precision check. *)
gminsol1=(<<solution2)[[1]]; {gminsol,gtmin}=(<<solution1)[[{1,2}]];

(* Check to the time beyond T upto which they are closer than tolerance 'tol' for each variable in the system. *)
(* This comparison should be made only till the time tempT upto which the candidate solution passed the accuracy check. *)
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."];];

(* Also display the actual tolerance achieved over specified time interval, T. *)
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<>"*"]];