diff --git a/Models20to30/LifeCycleModel21.m b/Models20to30/LifeCycleModel21.m
index f22c798e..ab4b80fd 100644
--- a/Models20to30/LifeCycleModel21.m
+++ b/Models20to30/LifeCycleModel21.m
@@ -77,7 +77,7 @@
 %% Grids
 % The ^3 means that there are more points near 0 and near 10. We know from theory that the value function will be more 'curved' near zero assets,
 % and putting more points near curvature (where the derivative changes the most) increases accuracy of results.
-a_grid=10*(linspace(0,1,n_a).^3)'; % The ^3 means most points are near zero, which is where the derivative of the value fn changes most.
+a_grid=15*(linspace(0,1,n_a).^3)'; % The ^3 means most points are near zero, which is where the derivative of the value fn changes most.
 
 % Grid for labour choice
 h_grid=linspace(0,1,n_d)'; % Notice that it is imposing the 0<=h<=1 condition implicitly
@@ -137,6 +137,9 @@
 end
 Params.mewj=Params.mewj./sum(Params.mewj); % Normalize to one
 AgeWeightsParamNames={'mewj'}; % So VFI Toolkit knows which parameter is the mass of agents of each age
+
+simoptions.optimize_nProbs=1;
+simoptions.verbose=2;
 StationaryDist=StationaryDist_FHorz_Case1(jequaloneDist,AgeWeightsParamNames,Policy,n_d,n_a,n_z,N_j,pi_z,Params,simoptions);
 % Note: Because we have simoptions.ExogShockFn, what we input for z_grid and  pi_z will just be ignored.
 
@@ -165,7 +168,7 @@
 title('Life Cycle Profile: Fraction Time Worked (h)')
 subplot(5,1,2); plot(1:1:Params.J,AgeConditionalStats.earnings.Mean)
 title('Life Cycle Profile: Labor Earnings (w kappa_j z h)')
-subplot(5,1,3); plot(1:1:Params.J,AgeConditionalStats.assets.Mean)
+subplot(5,1,3); hold on; plot(1:1:Params.J,AgeConditionalStats.assets.Mean); plot(1:1:Params.J,AgeConditionalStats.assets.Maximum); hold off
 title('Life Cycle Profile: Assets (a)')
 subplot(5,1,4); plot(1:1:Params.J,[AgeConditionalStats.fractionunemployed.Mean(1:Params.Jr-1),nan(1,Params.J-Params.Jr+1)])
 title('Life Cycle Profile: Fraction Unemployment (z==0)')
@@ -194,13 +197,18 @@
 ReturnFn=@(h,aprime,a,z,w,sigma,psi,eta,agej,Jr,pension,r,kappa_j,wg1,wg2,wg3,beta,sj) ...
     LifeCycleModel8_ReturnFn(h,aprime,a,z,w,sigma,psi,eta,agej,Jr,pension,r,kappa_j,wg1,wg2,wg3,beta,sj);
 [V_nomedical, Policy_nomedical]=ValueFnIter_Case1_FHorz(n_d,n_a,n_z,N_j, d_grid, a_grid, z_grid, pi_z, ReturnFn, Params, DiscountFactorParamNames, [], vfoptions);
+simoptions.optimize_nProbs=1;
+simoptions.verbose=2;
 StationaryDist_nomedical=StationaryDist_FHorz_Case1(jequaloneDist,AgeWeightsParamNames,Policy_nomedical,n_d,n_a,n_z,N_j,pi_z,Params,simoptions);
 AgeConditionalStats_nomedical=LifeCycleProfiles_FHorz_Case1(StationaryDist_nomedical,Policy_nomedical,FnsToEvaluate,Params,[],n_d,n_a,n_z,N_j,d_grid,a_grid,z_grid,simoptions);
 
 figure(2)
+hold on
 plot(1:1:Params.J,AgeConditionalStats.assets.Mean,1:1:Params.J,AgeConditionalStats_nomedical.assets.Mean)
+plot(1:1:Params.J,AgeConditionalStats.assets.Maximum,1:1:Params.J,AgeConditionalStats_nomedical.assets.Maximum)
+hold off
 title('Life Cycle Profile: Assets (a)')
-legend('Medical Expense Shocks','No Medical Shocks')
+legend('Medical Expense Shocks','No Medical Shocks','Medical Expense Shocks (Maximum)','No Medical (Maximum)')
 % Notice that medical expense shocks late in life cause elderly households
 % to hold more assets (as self-insurance against medical expense shocks)
 
diff --git a/Models20to30/LifeCycleModel23.m b/Models20to30/LifeCycleModel23.m
index d537505b..af5ba0bd 100644
--- a/Models20to30/LifeCycleModel23.m
+++ b/Models20to30/LifeCycleModel23.m
@@ -206,6 +206,10 @@
 end
 Params.mewj=Params.mewj./sum(Params.mewj); % Normalize to one
 AgeWeightsParamNames={'mewj'}; % So VFI Toolkit knows which parameter is the mass of agents of each age
+
+simoptions.optimize_nProbs=1;
+simoptions.verbose=2;
+simoptions.tolerance=1e-7;
 StationaryDist=StationaryDist_Case1_FHorz_PType(jequaloneDist,AgeWeightsParamNames,PTypeDistParamNames,Policy,n_d,n_a,n_z,N_j,N_i,pi_z,Params,simoptions);
 
 %% FnsToEvaluate are how we say what we want to graph the life-cycles of
diff --git a/Models20to30/LifeCycleModel28.m b/Models20to30/LifeCycleModel28.m
index 7cb2dcaf..e263eb5e 100644
--- a/Models20to30/LifeCycleModel28.m
+++ b/Models20to30/LifeCycleModel28.m
@@ -79,6 +79,7 @@
 simoptions.gridinterplayer=vfoptions.gridinterplayer; % grid interpolation layer must also be set in simoptions (because it changes Policy size/interpretation)
 simoptions.ngridinterp=vfoptions.ngridinterp;
 
+
 % At the bottom of this code/script there are some lines showing you what
 % pi_semiz_J which is created internally looks like.
 
@@ -131,7 +132,7 @@
 %% Grids
 % The ^3 means that there are more points near 0 and near 10. We know from theory that the value function will be more 'curved' near zero assets,
 % and putting more points near curvature (where the derivative changes the most) increases accuracy of results.
-a_grid=10*(linspace(0,1,n_a).^3)'; % The ^3 means most points are near zero, which is where the derivative of the value fn changes most.
+a_grid=12*(linspace(0,1,n_a).^3)'; % The ^3 means most points are near zero, which is where the derivative of the value fn changes most.
 
 % First, the AR(1) process z
 [z_grid,pi_z]=discretizeAR1_FarmerToda(0,Params.rho_z,Params.sigma_epsilon_z,n_z);
@@ -223,6 +224,8 @@
 % Because evaluating pi_semiz_J requires the d_grid we also have to provide
 simoptions.d_grid=d_grid;
 
+simoptions.optimize_nProbs=1;
+simoptions.verbose=2;
 StationaryDist=StationaryDist_FHorz_Case1(jequaloneDist,AgeWeightsParamNames,Policy,n_d,n_a,n_z,N_j,pi_z,Params,simoptions);
 
 %% FnsToEvaluate are how we say what we want to graph the life-cycles of
@@ -250,7 +253,7 @@
 title('Life Cycle Profile: Fraction Time Worked (h)')
 subplot(5,1,2); plot(1:1:Params.J,AgeConditionalStats.earnings.Mean)
 title('Life Cycle Profile: Labor Earnings (w kappa_j z h)')
-subplot(5,1,3); plot(1:1:Params.J,AgeConditionalStats.assets.Mean)
+subplot(5,1,3); hold on; plot(1:1:Params.J,AgeConditionalStats.assets.Mean); plot(1:1:Params.J,AgeConditionalStats.assets.Maximum); hold off
 title('Life Cycle Profile: Assets (a)')
 subplot(5,1,4); plot(1:1:Params.J,AgeConditionalStats.ninfants.Mean)
 title('Life Cycle Profile: Number of Infants (n1)')
diff --git a/Models31to35/LifeCycleModel31.m b/Models31to35/LifeCycleModel31.m
index cccbce32..ec24cf28 100644
--- a/Models31to35/LifeCycleModel31.m
+++ b/Models31to35/LifeCycleModel31.m
@@ -111,7 +111,7 @@
 %% Grids
 % The ^3 means that there are more points near 0 and near 10. We know from theory that the value function will be more 'curved' near zero assets,
 % and putting more points near curvature (where the derivative changes the most) increases accuracy of results.
-a_grid=13*(linspace(0,1,n_a).^3)'; % The ^3 means most points are near zero, which is where the derivative of the value fn changes most.
+a_grid=16*(linspace(0,1,n_a).^3)'; % The ^3 means most points are near zero, which is where the derivative of the value fn changes most.
 
 % First, the AR(1) process z
 [z_grid,pi_z]=discretizeAR1_FarmerToda(0,Params.rho_z,Params.sigma_epsilon_z,n_z);
@@ -266,6 +266,9 @@
 end
 Params.mewj=Params.mewj./sum(Params.mewj); % Normalize to one
 AgeWeightsParamNames={'mewj'}; % So VFI Toolkit knows which parameter is the mass of agents of each age
+
+simoptions.optimize_nProbs=1;
+simoptions.verbose=2;
 StationaryDist=StationaryDist_FHorz_Case1(jequaloneDist,AgeWeightsParamNames,Policy,n_d,n_a,n_z,N_j,pi_z,Params,simoptions);
 % riskyasset requires the grids when simulating the agent distribution to be able to handle aprime(d,u). The grids are passed in simoptions.
 
@@ -293,7 +296,7 @@
 title('Life Cycle Profile: Share of savings invested in risky asset (riskyshare)')
 subplot(3,1,2); plot(1:1:Params.J,AgeConditionalStats.earnings.Mean)
 title('Life Cycle Profile: Labor Earnings (w kappa_j z)')
-subplot(3,1,3); plot(1:1:Params.J,AgeConditionalStats.assets.Mean)
+subplot(3,1,3); hold on; plot(1:1:Params.J,AgeConditionalStats.assets.Mean); plot(1:1:Params.J,AgeConditionalStats.assets.Maximum); hold off
 title('Life Cycle Profile: Assets (a)')
 
 
diff --git a/Models31to35/LifeCycleModel35.m b/Models31to35/LifeCycleModel35.m
index 12903fe1..0cfabf74 100644
--- a/Models31to35/LifeCycleModel35.m
+++ b/Models31to35/LifeCycleModel35.m
@@ -44,7 +44,6 @@
 simoptions.riskyasset=1;
 % When there is more than one endogenous state, the riskyasset is the last one
 
-
 % Specify Epstein-Zin preferences
 vfoptions.exoticpreferences='EpsteinZin';
 vfoptions.EZpositiveutility=0; % Epstein-Zin preferences in utility-units have to be handled differently depending on whether the utility function is positive or negative valued (this is all done internally, you just need to use vfoptions to specify which)
@@ -122,8 +121,8 @@
 %% Grids
 % The ^3 means that there are more points near 0 and near 10. We know from theory that the value function will be more 'curved' near zero assets,
 % and putting more points near curvature (where the derivative changes the most) increases accuracy of results.
-asset_grid=-3+13*(linspace(0,1,n_a(2)))'; % Note, I use equal spacing (normally would put most points near zero)
-% note: will go from -3 to 13-3
+asset_grid=-3+23*(linspace(0,1,n_a(2)))'; % Note, I use equal spacing (normally would put most points near zero)
+% note: will go from -3 to 23-3
 % Make it so that there is a zero assets
 % Find closest to zero assets
 [~,zeroassetindex]=min(abs(asset_grid));
@@ -294,6 +293,9 @@
 end
 Params.mewj=Params.mewj./sum(Params.mewj); % Normalize to one
 AgeWeightsParamNames={'mewj'}; % So VFI Toolkit knows which parameter is the mass of agents of each age
+
+simoptions.optimize_nProbs=1;
+simoptions.verbose=2;
 StationaryDist=StationaryDist_FHorz_Case1(jequaloneDist,AgeWeightsParamNames,Policy,n_d,n_a,n_z,N_j,pi_z,Params,simoptions);
 % riskyasset requires the grids when simulating the agent distribution to be able to handle aprime(d,u). The grids are passed in simoptions.
 
@@ -322,7 +324,7 @@
 title('Life Cycle Profile: Share of savings invested in risky asset (riskyshare)')
 subplot(4,1,2); plot(1:1:Params.J,AgeConditionalStats.earnings.Mean)
 title('Life Cycle Profile: Labor Earnings (w kappa_j z)')
-subplot(4,1,3); plot(1:1:Params.J,AgeConditionalStats.assets.Mean)
+subplot(4,1,3); hold on; plot(1:1:Params.J,AgeConditionalStats.assets.Mean); plot(1:1:Params.J,AgeConditionalStats.assets.Maximum); hold off
 title('Life Cycle Profile: Assets (a)')
 subplot(4,1,4); plot(1:1:Params.J,AgeConditionalStats.housing.Mean)
 title('Life Cycle Profile: Housing (h)')
diff --git a/Models40to44/LifeCycleModel41.m b/Models40to44/LifeCycleModel41.m
index 7c730e74..ca300f08 100644
--- a/Models40to44/LifeCycleModel41.m
+++ b/Models40to44/LifeCycleModel41.m
@@ -94,9 +94,9 @@
 %% Grids
 % The ^3 means that there are more points near 0 and near 10. We know from theory that the value function will be more 'curved' near zero assets,
 % and putting more points near curvature (where the derivative changes the most) increases accuracy of results.
-asset_grid=10*(linspace(0,1,n_a(1)).^3)'; % The ^3 means most points are near zero, which is where the derivative of the value fn changes most.
+asset_grid=30*(linspace(0,1,n_a(1)).^3)'; % The ^3 means most points are near zero, which is where the derivative of the value fn changes most.
 
-h_grid=linspace(0.1,3,n_a(2))'; % Because h is an experienceasset, it will be interpolated onto this grid and so we need less grid points than usual
+h_grid=linspace(0.1,2,n_a(2))'; % Because h is an experienceasset, it will be interpolated onto this grid and so we need less grid points than usual
 % Note: deliberately omit 0 from h_grid
 
 % First, the AR(1) process z
@@ -225,6 +225,9 @@
 end
 Params.mewj=Params.mewj./sum(Params.mewj); % Normalize to one
 AgeWeightsParamNames={'mewj'}; % So VFI Toolkit knows which parameter is the mass of agents of each age
+
+simoptions.optimize_nProbs=1;
+simoptions.verbose=2;
 StationaryDist=StationaryDist_FHorz_Case1(jequaloneDist,AgeWeightsParamNames,Policy,n_d,n_a,n_z,N_j,pi_z,Params,simoptions);
 
 %% FnsToEvaluate are how we say what we want to graph the life-cycles of
@@ -250,9 +253,9 @@
 title('Life Cycle Profile: Female Earnings (w h z p)')
 subplot(5,1,3); plot(Params.agejshifter+(1:1:Params.J),AgeConditionalStats.maleearnings.Mean)
 title('Life Cycle Profile: Male Earnings (y_m)')
-subplot(5,1,4); plot(Params.agejshifter+(1:1:Params.J),AgeConditionalStats.assets.Mean)
+subplot(5,1,4); hold on; plot(Params.agejshifter+(1:1:Params.J),AgeConditionalStats.assets.Mean); plot(Params.agejshifter+(1:1:Params.J),AgeConditionalStats.assets.Maximum); hold off
 title('Life Cycle Profile: Assets (a)')
-subplot(5,1,5); plot(Params.agejshifter+(1:1:Params.J),AgeConditionalStats.femaleLFPH.Mean)
+subplot(5,1,5); hold on; plot(Params.agejshifter+(1:1:Params.J),AgeConditionalStats.femaleLFPH.Mean); plot(Params.agejshifter+(1:1:Params.J),AgeConditionalStats.femaleLFPH.Maximum); hold off
 title('Life Cycle Profile: female labor force participation history (h)')
 % Remember, we added a huge childcare cost, so women choose not to work
 % during ages 27-31. You can see that their human capital (femaleLFPH) has
diff --git a/Models40to44/LifeCycleModel42.m b/Models40to44/LifeCycleModel42.m
index 76d65185..16e7d775 100644
--- a/Models40to44/LifeCycleModel42.m
+++ b/Models40to44/LifeCycleModel42.m
@@ -64,9 +64,9 @@
 
 % Uncertain Human Capital
 Params.alpha_h=0.6; % returns to scale of human capital production fn
-Params.delta_h=0.02; % depreciation rate of human capital
-Params.ability=3; % ability to produce human capital
-Params.sscaler=5; % scale from units of study-time into something more appropriate for the human capital production function
+Params.delta_h_j=0.02*[ones(1,31), linspace(1,2,46-31), linspace(2,4,61-46), 4*ones(1,81-61)]; % depreciation rate of human capital
+Params.ability_j=[linspace(3,1,21), ones(1,81-21)]; % ability to produce human capital
+Params.sscaler=0.5; % scale from units of study-time into something more appropriate for the human capital production function
 % i.i.d. shocks that make human capital uncertain
 Params.sigma_epsilon_u=0.01;
 
@@ -94,13 +94,17 @@
 Params.wg2=3; % degree to which bequests are a luxury good (>=1; =1 would be a normal good)
 Params.wg3=Params.sigma; % By using the same curvature as the utility of consumption it makes it much easier to guess appropriate parameter values for the warm glow
 
+% Tuition
+% When tuition is more expensive (8), agents get advanced degree later, when they can afford it
+% When tutition is less expensive (6), they hurry up and study while kappa_j is low
+Params.tuition=7;
 
 %% Grids
 % The ^3 means that there are more points near 0 and near 10. We know from theory that the value function will be more 'curved' near zero assets,
 % and putting more points near curvature (where the derivative changes the most) increases accuracy of results.
-asset_grid=10*(linspace(0,1,n_a(1)).^3)'; % The ^3 means most points are near zero, which is where the derivative of the value fn changes most.
+asset_grid=30*(linspace(0,1,n_a(1)).^3)'; % The ^3 means most points are near zero, which is where the derivative of the value fn changes most.
 
-h_grid=linspace(1,4,n_a(2))'; % Because h is an experienceasset, it will be interpolated onto this grid and so we need less grid points than usual
+h_grid=linspace(1,20,n_a(2))'; % Because h is an experienceasset, it will be interpolated onto this grid and so we need less grid points than usual
 % Note: deliberately omit 0 from h_grid
 
 % First, the AR(1) process z
@@ -130,7 +134,7 @@
 % the experience asset).
 
 % aprimeFn gives the value of hprime
-vfoptions.aprimeFn=@(s,h,u,alpha_h, delta_h, ability,sscaler) u*(ability*(h*(sscaler*s))^alpha_h+h*(1-delta_h));
+vfoptions.aprimeFn=@(s,h,u,alpha_h, delta_h_j, ability_j,sscaler) u*(ability_j*(h*(sscaler*s))^alpha_h+h*(1-delta_h_j));
 % The first three inputs must be (d,a,u) [in the sense of aprime(d,a,u)], then any parameters
 
 % We also need to tell simoptions about the experienceassetu
@@ -162,12 +166,12 @@
 simoptions.gridinterplayer=vfoptions.gridinterplayer; % grid interpolation layer must also be set in simoptions (because it changes Policy size/interpretation)
 simoptions.ngridinterp=vfoptions.ngridinterp;
 
-% To better understand the human capital production function, here is a graph of it
+% To better understand the human capital production function, here is a graph of it (at age 1)
 hprime_shu=zeros([n_a,n_u]);
 for hh=1:n_a(2)
     for ss=1:n_d(2)
         for uu=1:n_u
-            hprime_shu(ss,hh,uu)=vfoptions.aprimeFn(s_grid(ss),h_grid(hh),u_grid(uu),Params.alpha_h,Params.delta_h,Params.ability,Params.sscaler);
+            hprime_shu(ss,hh,uu)=vfoptions.aprimeFn(s_grid(ss),h_grid(hh),u_grid(uu),Params.alpha_h,Params.delta_h_j(1),Params.ability_j(1),Params.sscaler);
         end
     end
 end
@@ -187,8 +191,8 @@
 DiscountFactorParamNames={'beta','sj'};
 
 % Use 'LifeCycleModel42_ReturnFn'
-ReturnFn=@(l,s,aprime,a,h,z,w,sigma,eta,psi,agej,Jr,pension,r,wg1,wg2,wg3,beta,sj)...
-    LifeCycleModel42_ReturnFn(l,s,aprime,a,h,z,w,sigma,eta,psi,agej,Jr,pension,r,wg1,wg2,wg3,beta,sj);
+ReturnFn=@(l,s,aprime,a,h,z,w,sigma,eta,psi,agej,Jr,pension,r,wg1,wg2,wg3,beta,sj,tuition)...
+    LifeCycleModel42_ReturnFn(l,s,aprime,a,h,z,w,sigma,eta,psi,agej,Jr,pension,r,wg1,wg2,wg3,beta,sj,tuition);
 % Notice how we have (l,s,aprime,a,h,z,...)
 % Follow same decision-next endo-endo-exo ordering as usual, but because h
 % is an experienceassetu, we do not include hprime as it is not chosen
@@ -257,7 +261,7 @@
 %% Initial distribution of agents at birth (j=1)
 % Before we plot the life-cycle profiles we have to define how agents are at age j=1. We will give them all zero assets.
 jequaloneDist=zeros([n_a,n_z],'gpuArray'); % Put no households anywhere on grid
-jequaloneDist(1,4,floor((n_z+1)/2))=1; % All agents start with zero assets, h_grid(4) of human capital, and the median shock [h_grid(4) is roughly same as y_m(1), just my arbitrary decision]
+jequaloneDist(1,2,floor((n_z+1)/2))=1; % All agents start with zero assets, h_grid(2) of human capital, and the median shock [h_grid(2) is roughly same as y_m(1), just my arbitrary decision]
 
 %% We now compute the 'stationary distribution' of households
 % Start with a mass of one at initial age, use the conditional survival
@@ -269,6 +273,9 @@
 end
 Params.mewj=Params.mewj./sum(Params.mewj); % Normalize to one
 AgeWeightsParamNames={'mewj'}; % So VFI Toolkit knows which parameter is the mass of agents of each age
+
+simoptions.optimize_nProbs=1;
+simoptions.verbose=2;
 StationaryDist=StationaryDist_FHorz_Case1(jequaloneDist,AgeWeightsParamNames,Policy,n_d,n_a,n_z,N_j,pi_z,Params,simoptions);
 
 
@@ -295,9 +302,9 @@
 title('Life Cycle Profile: Study time (s)')
 subplot(5,1,3); plot(Params.agejshifter+(1:1:Params.J),AgeConditionalStats.earnings.Mean)
 title('Life Cycle Profile: Earnings (w h z l)')
-subplot(5,1,4); plot(Params.agejshifter+(1:1:Params.J),AgeConditionalStats.humancapital.Mean)
+subplot(5,1,4); hold on; plot(Params.agejshifter+(1:1:Params.J),AgeConditionalStats.humancapital.Mean); plot(Params.agejshifter+(1:1:Params.J),AgeConditionalStats.humancapital.Maximum); hold off
 title('Life Cycle Profile: Human Capital (h)')
-subplot(5,1,5); plot(Params.agejshifter+(1:1:Params.J),AgeConditionalStats.assets.Mean)
+subplot(5,1,5); hold on; plot(Params.agejshifter+(1:1:Params.J),AgeConditionalStats.assets.Mean); plot(Params.agejshifter+(1:1:Params.J),AgeConditionalStats.assets.Maximum); hold off
 title('Life Cycle Profile: Assets (a)')
 % Because the calibration is a bit silly, assets are an inferior investment
 % to human capital. Hence throughout working life households just invest in
diff --git a/Models40to44/LifeCycleModel42_ReturnFn.m b/Models40to44/LifeCycleModel42_ReturnFn.m
index 61fcea76..13cfbb2e 100644
--- a/Models40to44/LifeCycleModel42_ReturnFn.m
+++ b/Models40to44/LifeCycleModel42_ReturnFn.m
@@ -1,12 +1,12 @@
-function F=LifeCycleModel42_ReturnFn(l,s,aprime,a,h,z,w,sigma,eta,psi,agej,Jr,pension,r,wg1,wg2,wg3,beta,sj)
+function F=LifeCycleModel42_ReturnFn(l,s,aprime,a,h,z,w,sigma,eta,psi,agej,Jr,pension,r,wg1,wg2,wg3,beta,sj,tuition)
 
 leisure=1-s-l; % one unit of time, minus time spent studying and working
 
 F=-Inf;
 if agej<Jr % If working age
-    c=w*h*z*l+(1+r)*a-aprime; % Add z here
+    c=w*h*z*l+(1+r)*a-aprime-s*tuition; % Add z here
 else % Retirement
-    c=pension+(1+r)*a-aprime;
+    c=pension+(1+r)*a-aprime-s*tuition/2;
 end
 
 if c>0 && leisure>0
diff --git a/ModelsAppendixA/LifeCycleModelA6.m b/ModelsAppendixA/LifeCycleModelA6.m
index ac3cbb71..deb527bd 100644
--- a/ModelsAppendixA/LifeCycleModelA6.m
+++ b/ModelsAppendixA/LifeCycleModelA6.m
@@ -184,6 +184,7 @@
 simoptions.gridinterplayer=vfoptions.gridinterplayer; % grid interpolation layer must also be set in simoptions (because it changes Policy size/interpretation)
 simoptions.ngridinterp=vfoptions.ngridinterp;
 
+
 %% Now, create the return function
 DiscountFactorParamNames={'beta','sj'};
 
@@ -234,6 +235,9 @@
 end
 Params.mewj=Params.mewj./sum(Params.mewj); % Normalize to one
 AgeWeightsParamNames={'mewj'}; % So VFI Toolkit knows which parameter is the mass of agents of each age
+
+simoptions.optimize_nProbs=1;
+simoptions.verbose=2;
 StationaryDist=StationaryDist_FHorz_Case1(jequaloneDist,AgeWeightsParamNames,Policy,n_d,n_a,n_z,N_j,pi_z,Params,simoptions);
 
 
@@ -260,7 +264,7 @@
 title('Life Cycle Profile: Fraction Time Worked (h)')
 subplot(3,1,2); plot(1:1:Params.J,AgeConditionalStats.earnings.Mean)
 title('Life Cycle Profile: Labor Earnings (w kappa_j h e1 e2)')
-subplot(3,1,3); plot(1:1:Params.J,AgeConditionalStats.assets.Mean)
+subplot(3,1,3); hold on; plot(1:1:Params.J,AgeConditionalStats.assets.Mean); subplot(3,1,3); plot(1:1:Params.J,AgeConditionalStats.assets.Maximum); plot(1:1:Params.J,AgeConditionalStats.assets.Minimum); hold off
 title('Life Cycle Profile: Assets (a)')