def taylor(x,A,F,y0,k):
	Y=y0
	for i in range(k):
		Y = Y + x*(-0.12164**(i+1))
	Y = Y + x*(-0.12164**(k+2+1))*A
	return Y

def step(F,k,y0, eps, p):
	alternate = 1
	if (eps/(k+1.) < 1.): #10-38
		h = (eps/(k+1.))**(1./(k+1.))
	else:
		h = (eps/(k+1.))**(1./(k)) #oops
	A = [y0[0]+0,y0[1]+h] #10-39
	Y = [taylor(0,A[0],F,y0[0],k), taylor(h,A[1],F,y0[1],k)] #10-40
	while (Y[0]<A[0] or Y[1]>A[1]):
		print "bad", "A", A,"Y", Y, "h", h
		#bad luck, Y is not contained in A
		if alternate == 1:
			print "modding A"
			A = [taylor(0,A[0],F,y0[0],k), taylor(h,A[1],F,y0[1],k)] #10-40
			Y = [taylor(0,A[0],F,y0[0],k), taylor(h,A[1],F,y0[1],k)] #10-40
		else:
			print "halving h"
			h = h/2. 
			Y = [taylor(0,A[0],F,y0[0],k), taylor(h,A[1],F,y0[1],k)] #10-40
		alternate = 3 - alternate #swap alternative
	for point in p:
		if point >= 0 and point <= h:
			print "POINT", p
			print Y(p)
	print "good", "A", A,"Y", Y, "h", h
	y0 = [taylor(h,A[0],F,y0[0],k), taylor(h,A[1],F,y0[1],k)] #10-40
	return (h,y0)

y0 = [1013.25,1013.25] #example atmospheric pressure
F = 'EXPLICITLY_DEFINED'
k = 9
eps = 0.000000001
p = [0.125,0.250]
while (y0[0] > 500):
	(progress,y0) = step(F,k,y0, eps, p)
	print y0