Below is what I was able to come up with. I checked the first couple iterations to make sure they were correct. After that, it looks right at a glance, but I haven't verified. This is from my old SAGE 8.4 installation (I have 9.0 on a different computer), so the `print`

statement will be formatted differently in version 9.0.

The key is to substitute the result back into `v0`

*before* substituting back into `vk`

.

```
sage: v0 = log(k)
sage: for n in range(5):
....: vk = log(k-k1) + beta*v0(k=k1)
....: FOC = vk.diff(k1)
....: k1star = solve(FOC==0, k1)
....: print(n, k1star)
....: v0 = (vk).subs(k1=k1star[0].rhs())
....:
(0, [k1 == beta*k/(beta + 1)])
(1, [k1 == (beta^2 + beta)*k/(beta^2 + beta + 1)])
(2, [k1 == (beta^3 + beta^2 + beta)*k/(beta^3 + beta^2 + beta + 1)])
(3, [k1 == (beta^4 + beta^3 + beta^2 + beta)*k/(beta^4 + beta^3 + beta^2 + beta + 1)])
(4, [k1 == (beta^5 + beta^4 + beta^3 + beta^2 + beta)*k/(beta^5 + beta^4 + beta^3 + beta^2 + beta + 1)])
```

I don't see the point : the value of

`k1`

that nullifies $\frac{\partial vk}{\partial k1}$ isby definitionan expression that doesnotdepend on`k1`

.At this point :

But, of course :

So there's

nothingto iterate.Aren't you trying to implement something like Newton's method of root-finding when

`solve`

doesn't give a explicit answer ?The context is given in Ask Sage question 52559: Cake-eating problem.

@EmmanuelCharpentier I see the issue with my approach. Though I actually made a video on solving a planners problem which is similar to this that uses value function iteration, https://www.youtube.com/watch?v=Bv_tP.... I just want to be able to code something like this.