Thank you. Well, I think that in the example the production function is y = k^\alpha, but the analytical functions are for (k-c)^\alpha, i.e. for the problem that is solved.
From the book, we have:

So: the problem tweaks are as follows:
In the class OptimalGrowthModel
:
def state_action_value(self, c, y, v_array):
"""
Right hand side of the Bellman equation.
"""
u, f, β = self.u, self.f, self.β
v = interp1d(self.grid, v_array, fill_value="extrapolate") # correction -> fill_value="extrapolate"
return u(c) + β * v( f(y) - c ) # correction -> v(f(y) - c)
Then we need to change analytical function v_star
, to say v_super_star
:
def v_super_star(k, α, β, μ):
"""
True value function for A * (k**α) case
"""
A = 1
c1 = np.log(A * (1 - α * β) )
c2 = α * β * np.log(A * α * β)
c3 = 1/(1 - β)
c4 = 1/(1 - α * β)
return c3 * (c1 + c2 * c4) + c4 * α * np.log(k)
Finally, add another plot to see that the values have changed.
grid = og.grid
v_init = v_star(grid, α, og.β, og.μ) # Start at the solution
v_init_1 = v_super_star(grid, α, og.β, og.μ) # correction
v_greedy, v = T(v_init_1, og) # Apply T once to correction
fig, ax = plt.subplots()
ax.set_ylim(-35, -24)
ax.plot(grid, v, lw=2, alpha=0.6, label='$Tv_{update}^*$') # name correction
ax.plot(grid, v_init, lw=2, alpha=0.6, label='$v^*$')
ax.plot(grid, v_init_1, lw=2, alpha=0.6, label='$v_{update}^*$') # correction
ax.legend()
plt.show()