function splineList = FindDPPolicy(utilFun, gamma, wealthGrid, ...
    incomeValues, transMatrix, retScenarios, retProbs, riskFree)
timeHorizon = length(wealthGrid);
numScenarios = length(retProbs);
numEmplStates = length(incomeValues);
% splines will be contained in a cell array
splineList = cell(timeHorizon, numEmplStates);
% last step, consume all available wealth, plus last income
for j = 1:numEmplStates
    wealthValues = wealthGrid{timeHorizon};
    auxFunValues = feval(utilFun,wealthValues+incomeValues(j));
    splineList{timeHorizon, j}= spline(wealthValues, auxFunValues);
end  % for employment states
% Now carry out backward recursion
options = optimoptions(@fmincon,'Display','off');
for t = (timeHorizon-1):-1:1
    fprintf(1, 'Working on time period %d\n',t);
    wealthValues = wealthGrid{timeHorizon};
    numPoints = length(wealthValues);
    for j = 1:numEmplStates
        for i = 1:numPoints
            x0 = [wealthValues(i); 0.5]; % starting point for fmincon
            [~,outVal] = fmincon(@(x) -objfun(x,i,j),x0,[],[],[],[], ...
                zeros(2,1), [wealthValues(i)+incomeValues(j);1], ...
                [],options);
            auxFunValues(i) = -outVal;
        end % for wealth gridpoints
        % build spline for this employment state and time
        splineList{t,j}= spline(wealthValues, auxFunValues);
    end % for employment states
end %for time
% nested subfunction
    function outVal = objfun(x,idxW,idxE)
        C = x(1);
        alpha = x(2);
        outVal = 0; % expected value function for future states
        for hh=1:numScenarios
            futureFinWealth = (wealthValues(idxW)+incomeValues(idxE)-C)* ...
                (1+riskFree+alpha*(retScenarios(hh)-riskFree));
            for jj=1:numEmplStates
                outVal = outVal + retProbs(hh)*transMatrix(idxE,jj) * ...
                    ppval(splineList{t+1,jj}, futureFinWealth);
            end % for jj
        end % for hh
        outVal = feval(utilFun,C)+gamma*outVal;
    end
end


