import swiginac
import sys

if len(sys.argv) < 4:
    print "Usage: \"python kthorder.py \'$ODE\' IC k\" where $ODE is the ODE with y as dependent variable and x as independent variable, IC the initial condition, and k the order of approximation"
    sys.exit(1)

x = swiginac.symbol('x')
y = swiginac.symbol('y')
F = eval(sys.argv[1])
IC = float(sys.argv[2])
k = int(sys.argv[3])

#does not work yet, unclear how to eval HOF in general

def taylor(xx,A,F,y0,k):
	Y = y0
	Yprime = F
	for i in range(k):
		Yprime = F.diff()
	print "taylor", F.series(x==0.0,k)
	return F.series(x==0.0, k).subs(x==A)

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 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 = [IC,IC] 
k = 9
eps = 0.000000001
p = [1.0,1.1,1.2,1.3,1.4]
while (y0[0] < 3.):
	(progress,y0) = step(F,k,y0, eps, p)
	print y0