from scipy.integrate import odeint
from scipy import arange
from numarray import *
from array import array
from math import *


tauvalues=[3,5,7,8,9,10]

T=[]
base=1.01
x=log(1e-15)/log(base)
T.append(base**x)
while T[len(T)-1] < 1e10 :
	x+=1.0
	T.append(base**x)
t=zeros(len(T)+1,"Float64")
for i in range(len(T)-1):
	t[i+1]=T[i]
last=len(t)-1
t[last]=1e10

k1=10**-9
k2=10**9
g1=g3=1
g2=(k1 / (k1-k2))

def f(y,x):
	return [-k1*y[0],k1*y[0] - k2*y[1],k2*y[1]]

f0=[1,0,0]

def y(t):
	return [g1 * exp(-k1*t),(g1 * exp(-k1*t) * (k1 / (k2-k1))) + (g2*exp(-k2*t)),(g1 * exp(-k1*t) * (k1 / (k1-k2))) - (g2*exp(-k2*t)) + g3]

def answers(sol):
	nfe=sol[1]['nfe']
	numeval=nfe[len(nfe)-1]
	nje=sol[1]['nje']
	numjeval=nje[len(nje)-1]
	#get answer at endpoint
	b=sol[0].tolist()
	b=b[len(b)-1]
	#give it back
	return b,numeval,numjeval

print "\\begin{tabular}{|lll|} \\hline"

for i in tauvalues:
	tau=10**-i
	#soln=odeint(f,f0,t,atol=0.001,rtol=tau,full_output=1)
	soln=odeint(f,f0,t,rtol=tau,full_output=1)
	soln=answers(soln)
	#print soln
	yn=soln[0]
	err=zeros(3,"Float64")
	for j in range(3):
		err[j] = abs(y(1e10)[j] - yn[j]) / abs(y(1e10)[j])

	#print output stuff
	print "$\\tau$ = ",tau," & func eval =",soln[1]," & jac evals = ",soln[2],"\\\\ \\hline"
	print "$y_{n}$ & rel err & \\\\"
	for j in range(3):
		print yn[j]," & ",err[j],"& \\\\"
	print "\\hline \\hline"

print "\\end{tabular}"