import itertools
import numpy as np
from scipy import optimize
from scipy import integrate

n=763.0; s=762.0; b=0.0026; y=0.5673
def p(x,b,n,y,s): return 1.0/(x*(b*n+y*np.log(x)-s*b*x))
def Ui(x,x1,b,n,y,s): return (integrate.fixed_quad(p, x, x1, args=(b,n,y,s),    n=5)[0]-1.0)

def t(x,x1,b,n,y,s): return optimize.newton(Ui, x1, p, args=(x,b,n,y,s),  maxiter=500)

x=0.99;  x1=1.0
X=[1]
for _ in xrange(1,15): x1=round(t(x1,x, b, n, y, s),500); x=x1; X.append(x)
Xu = X[1:]