% Session p. 128

f=inline('sum(OccupationProbability(epsilon,K))-2','epsilon','K');
for K=125*2.^(0:3)
    out=sprintf('K = %4i \t\t epsilon* = %17.15f',K,...
        fzero(f,[0.06,0.07],optimset('TolX',1e-14),K));
    disp(out);
end