This file implements the algorithm 'Forward Shoot 1D'  outlined in paper titled, "Smart Forward Shooting" which solves models with one jump variable 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 1 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),σ=(50/100),ρ=(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),ω=(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 *)
{
cdot=τ c (r-ρ);
kdot = Y[k]-c-δ k;
system  ={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]-δ;
};

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

(* INITIALIZATION CELL - Initialization for solving for the global saddle path. *)
{  
T = 40;                        (* solution time (years). BEWARE!!! Setting T larger than necessary can substantially increase computational time. *); 
tol = (10^(-4));                (* desired accurcy of the solutions of the variables in the system *);
precision = 16;     (* starting value of precision of calculations *) ;
epsilon= ((10^(-0))/40)(Subscript[c, o]); (* radius of the ϵ-segment *)
pini={(Subscript[c, o])};              (* center of the initial ϵ-setment *);
};
(* 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},{c,k},{ (Subscript[c, o]),(Subscript[k, o])},{(Subscript[c, f]),(Subscript[k, f])},{{r,rfunc,{k,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];
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,rfunc[(Subscript[k, o]),0],rfunc[k[t],t]]/.gminsol,rfunc[(Subscript[k, f]),Infinity]},{t,-tplot/10,tplot},PlotRange->{.099,.107},AxesOrigin->{-tplot/10,.099},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];

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

(* Tabulate the final solutions. *)
Labels={"c","i","k","Y","r"};InitialSteadyState={(Subscript[c, o]),δ (Subscript[k, o]),(Subscript[k, o]),(Subscript[Y, o]),ρ}//N;FinalSteadyState={(Subscript[c, f]),δ (Subscript[k, f]),(Subscript[k, f]),(Subscript[Y, f]),ρ}//N;
(Dat=Table[{tt,c[t],Y[k[t]]-c[t],k[t],Y[k[t]],rfunc[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 1D 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 variable 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 ϵ-segment to epsilon. Note that setting T larger than necessary can significantly increase computation time. \n
pini  = {a} gives center of the initial segment in the jump space. The typical value is {(Subscript[x, o])}, the initial steady state values of the jump variable 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,epsilon,Tmax,firstsol,precisionini,system,boxvars,box,pmin,digits,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 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:=Evaluate[Apply[Or,Join[MapThread[First[Map[(#1),{t}]]<((Min[#2,#3])/20)&,boxvars],MapThread[First[Map[(#1),{t}]]>5Max[#2,#3]&,boxvars]]]];

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

Label[begin];
$MinPrecision=digits=precision;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}=compdis[pmin,0];

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

(* SOLVE - Solve for the global nonlinear saddle path *);
{gmindis,gtmin,gminsol}=compdis[pmin,epsilon];

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

If[gtmin==(10^100),
	If[precision≥80,
		Print["The system is highly ill-conditioned."];Abort[],
		precision=precision+4;Print["Increased 'precision' by 4 digits to ",precision," digits."];Goto[begin]];];

];
 
(* 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_]:=Block[{dis,dislast,tdis,steps,tbeg,tend,intsz,traj,norm},
dis= (10^7);dislast=(10^8);steps=Max[4Floor[tsol+1],450] ;tbeg=0;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 *)
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_,ϵ_]:=Block[{iflag,it,shrink,pexp,difftmin,notsatistmin,satistmin,diffprec,notsatisprec,satisprec,increaseprec,step2,step3,congratulations,done,pini,ini,dis,tdis,jdis,jtdis,dissol,jp,sol,pmin,epsilon},

iflag=0;it=1;shrink=False;{pmin,epsilon}={p,ϵ};
If[Length[sysexpand]==0,pexp={},pexp=First[sysexpand][Apply[Sequence,Last[sysexpand]]];];

(* 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.";

Label[begin];
If[ϵ≠0,Print["\nBEGINNING ITERATION ",it," WITH CENTER OF SEGMENT AT ",N[pmin,precision]]];

For[j=-1,j≤ 1  ,j=j+2,
	pini=pmin+{j epsilon};
	ini=Thread[Table[vars[[i]][0],{i,1,Length[vars]}]==Apply[Sequence,{pini,pexp,Drop[varso,1+Length[pexp]]},1],List];
	sol=solvesys[ini,0,Tmax ];{dis,tdis}=distance[sol,sol[[1,1,2,1,1,2]]];jdis[j]=dis;jtdis[j]=tdis;dissol[j]=sol;jp[j]=pini;
];
	
If[ϵ==0,gminsol=sol;{gmindis,gtmin} ={dis,tdis};pmin=pini;Goto[end] ];
If[Min[jdis[-1],jdis[1]]>gmindis,shrink=True];

If[jdis[-1]<gmindis,
	{gmindis,gtmin,gminsol,pmin}={jdis[-1],jtdis[-1],dissol[-1],jp[-1]},
	
	If[jdis[1]<gmindis,
		{gmindis,gtmin,gminsol,pmin}={jdis[1],jtdis[1],dissol[1],jp[1]},
		
		If[shrink==True,epsilon=epsilon/2,epsilon=epsilon/1]]];

Print["AT THE END OF ITERATION = ",it,"     TIME OF CLOSEST APPROACH = ",gtmin,"     MINIMUM DISTANCE = ",gmindis];

(* ERROR CHECK - Check for the precision of the solution over the time period of interest. *) 
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,
			{gminsol,gtmin,gmindis,pmin,epsilon}>>solution2;
			Print[""];diffprec;Print[notsatisprec];precision=precision+4;increaseprec,
			
			{gminsol,gtmin,gmindis,pmin,5epsilon}>>solution1;
			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. *);
];

If[(epsilon/(pmin[[1]]))<(10^(-precision)),Print["You may be losing precision in calculations."];gtmin=(10^100);Goto[end];];

it++;Goto[begin];

Label[end];
{gmindis,gtmin,gminsol}];

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[{gminsol1,j},
(* 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. *)
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<>"*"]];