from scipy import *
from scipy.optimize import bisect
from scipy import gplt
import math

#here we solve the bvp y'' = y using shooting methods.

#globals
H = 1.0/50.0

def initCond(eta):
	ic = zeros( 2, Float64 )
	ic[0] = 1.0
	ic[1] = eta
	return ic

def euler(ic,h,left,right):
	yn = [ic[0]]
	un = [ic[1]]
	x = arange(left,right + h,h)

	for i in range(1,len(x)) :
		yn.append( yn[i-1] + h*un[i-1] )
		un.append( un[i-1] + h*yn[i-1] )

	return [x,yn]

#our numerical solution
def F(eta) :
	return euler(initCond(eta),H,-1,1)[1][100] - 1

#the true solution
def y(x) :
	return (math.exp(x+1) + math.exp(1-x)) / (math.exp(2) + 1)

#use these lines to determine bracket for bisection, just eyeball it
#p = arange(-1.0,1.0,0.1)
#for i in p :
#	print i,"&",F(i),"\\\\ \hline"

#find good value for eta, from previous step
trueEta = bisect(F,-0.8,-0.7)

#numeric solution for this problem
ans = euler(initCond(trueEta),H,-1,1)

#save a few other solutions for the writeup
p = arange(-1.0,-0.5,0.1)
plotSoln = [ans]
for i in p :
	plotSoln.append( euler(initCond(i),H,-1,1) )

trueY = []
for i in ans[0] :
	trueY.append( y(i) )

X = []
Y = []
for i in range(len(ans[0])) :
	#X.append( [ans[0][i],ans[0][i]] )
	#Y.append( [ans[1][i],trueY[i]] )
	t_x = []
	t_y = []
	for j in range( len(plotSoln) ) :
		t_x.append( plotSoln[j][0][i] )
		t_y.append( plotSoln[j][1][i] )
	X.append( t_x )
	Y.append( t_y )

	

#graph of solution
#gplt.plot(ans[0],ans[1]) #just the numerical solution
gplt.plot(X,Y) #both numerical and true solutions

#uncomment the following line to save a copy of the plot to img.png
gplt.output("img.png",'png color')