#Nathan Fox
#Homework 16
#I give permission for this work to be posted online

#Read procedures from class
read(`C12.txt`):
read(`C16.txt`):

Help:=proc():
    print(` BestAdams(f,x,y,x1,D:={1,2,3,4,5,6,7,8,9,10},H:={1/10,1/100,1/1000}) `):
end:

##PROBLEM 1##

#BestAdams(f,x,y,x1): inputs an expression f in x and y such that 
#dsolve({diff(y(x),x)=f,y(0)=1},y(x)); 
#gives you a closed-form answer, and outputs the d in D (default
#{1,2,3,4,5,6,7,8,9,10}) and h in H (default {1/10, 1/100, 1/1000})
#such that 
#AAB(f,x,y,0,1,h,x1,FAB(d)) 
#is as close as possible to the true y(x1),
#and also returns the error.
#(return value is an expression sequence d, h, error)
BestAdams:=proc(f,x,y,x1,D:={1,2,3,4,5,6,7,8,9,10},H:={1/10,1/100,1/1000}) local dsv, d, h, champ, rec, ab:
    dsv:=DS(f, x, y, 0, 1):
    if dsv = FAIL then
        return FAIL:
    fi:
    dsv:=evalf(subs(x=x1, dsv)):
    rec:=infinity:
    for d in D do
        for h in H do
            ab:=AAB(f,x,y,0,1,h,x1,FAB(d))[2][-1]:
            if abs(ab-dsv) < rec then
                champ:=d,h:
                rec:=abs(ab-dsv):
            end:
        od:
    od:
    return champ,rec:
end:

#BestAdams(y,x,y,1) returns d=2, h=1/1000, error=.2489e-5
#BestAdams(x+y,x,y,5) took too long
#BestAdams(x+y,x,y,5,{1,2,3,4,5,6,7,8,9,10},{1/10,1/100}) returns d=3, h=1/100, error=.300470e-1
#BestAdams(x+3*y,x,y,10) took too long
#BestAdams(x+3*y,x,y,10,{1,2,3,4,5,6,7,8,9,10},{1/10,1/100}) returns d=3, h=1/100, error=.1397307e11 (in other words, they all suck)

##PROBLEM 2##

#GIRK(RK4, -15*y, x, y, 0, 1, 0.1, 1) gives 0.97e-13 for y(1)
#GIRK(GL, -15*y, x, y, 0, 1, 0.1, 1) gives 0.34e-6 for y(1)
#The correct answer is exp(-15)=0.31e-6