I noticed that to compute the distribution using the Tan method (which is the default) we do not use parallel computing on CPU, which makes sense. But nonetheless the code automatically starts a parallel pool, which slows things down without any apparent benefit. The problem is in these lines of code:

if ~isfield(simoptions,'ncores')
        try
            PoolDetails=gcp;
            simoptions.ncores=PoolDetails.NumWorkers;
        catch
            simoptions.ncores=1;
        end
    end

which are around line 55 in StationaryDist_FHorz_Case1. I think that there is no need to start a parallel pool on the cpu: with the Tan improvement we do sparse matrices on the cpu, and when using gpuArray we don't need a parallel pool either.

In line 31 of ValueFnIter/InfHorz/ValueFnIter_Case1_NoD_Par2_raw.m, one can use automatic expansion (better for memory and speed)

The term EVz is of size (na,1) and is automatically expanded to (na,na), so I would eliminate the multiplication by ones(1,na)

Hi and thank you for this package!
Can your toolkit be used to solve default models (w/ nested value functions)?
For example: (Matlab Duarte etal)

In StationaryDist_FHorz_Case1_Iteration_TwoProbs_raw, I would replace the line 27 of the function

Gammatranspose=sparse(Policy_aprimez(:,:,jj),II2,PolicyProbs(:,:,jj),N_a*N_z,N_a*N_z);

with

Gammatranspose=sparse(Policy_aprimez(:,1,jj),II1,PolicyProbs(:,1,jj),N_a*N_z,N_a*N_z)+...
        sparse(Policy_aprimez(:,2,jj),II1,PolicyProbs(:,2,jj),N_a*N_z,N_a*N_z);

On line 21, precompute
II1=(1:1:N_a*N_z)';

instead of

II2=[1:1:N_a*N_z; 1:1:N_a*N_z]';

I propose to replace

seedpointvec=max(cumsumInitialDistVec>rand(1));

with

[~,seedpointvec]=max(cumsumInitialDistVec>rand(1));

wherever it occurs in the simpanelindexes code

Hello Robert
I have some problems when trying to re-run RestucciaRogerson2008.m. Here are some error messages:
******* Error Message Begin ********
Incorrect use of fieldnames
The input parameter of type 'cell' is invalid. The input must be a struct or a Java or COM object.

Erros in HeteroAgentStationaryEqm_Case1_EntryExit_subfn (Line 136)
if ~isempty(fieldnames(FnsToEvaluate)) % Note that the entry/exit aggregates are treated seperately, so it is possible for this to be empty

Errors in
HeteroAgentStationaryEqm_Case1_EntryExit>@(p)HeteroAgentStationaryEqm_Case1_EntryExit_subfn(p,n_d,n_a,n_s,pi_s,d_grid,a_grid,s_grid,ReturnFn,FnsToEvaluate,GeneralEqmEqns,Parameters,DiscountFactorParamNames,ReturnFnParamNames,FnsToEvaluateParamNames,GeneralEqmEqnParamNames,GEPriceParamNames,EntryExitParamNames,heteroagentoptions,simoptions,vfoptions) (Line 102)
GeneralEqmConditionsFnOpt=@(p) HeteroAgentStationaryEqm_Case1_EntryExit_subfn(p, n_d, n_a, n_s, pi_s, d_grid, a_grid, s_grid, ReturnFn, FnsToEvaluate, GeneralEqmEqns, Parameters, DiscountFactorParamNames, ReturnFnParamNames, FnsToEvaluateParamNames, GeneralEqmEqnParamNames, GEPriceParamNames, EntryExitParamNames, heteroagentoptions, simoptions, vfoptions);

Erros in fminsearch (Line 201)
fv(:,1) = funfcn(x,varargin{:});

Errors in HeteroAgentStationaryEqm_Case1_EntryExit (Line 116)
[p_eqm_vec,GeneralEqmConditions]=fminsearch(GeneralEqmConditionsFnOpt,p0,minoptions);

Erros in HeteroAgentStationaryEqm_Case1 (Line 117)
[p_eqm,p_eqm_index, GeneralEqmConditions]=HeteroAgentStationaryEqm_Case1_EntryExit(n_d, n_a, n_z, n_p, pi_z, d_grid, a_grid, z_grid, ReturnFn, FnsToEvaluate, GeneralEqmEqns, Parameters, DiscountFactorParamNames, ReturnFnParamNames, FnsToEvaluateParamNames, GeneralEqmEqnParamNames, GEPriceParamNames, EntryExitParamNames, heteroagentoptions, simoptions, vfoptions);

Errors in RestucciaRogerson2008 (Line 302 )
[p_eqm,p_eqm_index, GeneralEqmCondition]=HeteroAgentStationaryEqm_Case1(0, n_a, n_z, n_p, pi_z, [], a_grid, z_grid, ReturnFn, FnsToEvaluate, GeneralEqmEqns, Params, DiscountFactorParamNames, ReturnFnParamNames, FnsToEvaluateParamNames, GeneralEqmEqnParamNames, GEPriceParamNames,heteroagentoptions, simoptions, vfoptions, EntryExitParamNames);

****** Error Message End**********

Thank you for your help!