<1]; (* Find the value of X where the constraint begins to bind *) (* ########################################################*) BindingPoint = (beta^(1/-rho) )OmegaInvPrimeTm1[0]; (* For each value of xGrid, calculate the maximum utility and optimal c *) xGrid = Union[{BindingPoint},xGrid]; (* ########################################################*) FOCInvConstrcValuesTm1 = Table[ tmpResults = FindRoot[ uprime[c] == beta (OmegaInvPrimeTm1[xGrid[[i]]-c])^-rho ,{c,{.001,xGrid[[i]]+.01}}]; (* The FindRoot command returns {control->optimum} The code below extracts the two elements *) (* ########################################################*) cUnconstr = (c /. tmpResults); If[cUnconstr > xGrid[[i]],xGrid[[i]],cUnconstr] (* ########################################################*) ,{i,Length[xGrid]}]; cTm1FOCInvConstr = Interpolation[Transpose[{xGrid,FOCInvConstrcValuesTm1}],InterpolationOrder->1]; PlotcTm1IntExpFOCInvConstr = Plot[{cTm1FOCInvConstr[x]},{x,xMin,xMax} ,DisplayFunction->Identity ]; PlotComparecTm1DE = Show[PlotcTm1IntExpFOCInv,PlotcTm1IntExpFOCInvConstr ,DisplayFunction->$DisplayFunction ,ImageSize->{72 6.,72 6./GoldenRatio} ,PlotLabel->"Comparison of cTm1FOCInv[x] with Constrained" ]; Display[ParentDirectoryString<>"PlotComparecTm1DE.EPS",PlotComparecTm1DE,"EPS"];