|   | 1 |  initial version  | 
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)])
 Copyright Sage, 2010. Some rights reserved under creative commons license. Content on this site is licensed under a Creative Commons Attribution Share Alike 3.0 license.
 
                
                Copyright Sage, 2010. Some rights reserved under creative commons license. Content on this site is licensed under a Creative Commons Attribution Share Alike 3.0 license.