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

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

Help:=proc():
    print(` FindRKP(s,d,p,x,y) , FindRK(s,d) , FindRKR(s,d) `):
end:

##PROBLEM 1##

#FindRKP(s,d,p,x,y): inputs positive integers s and d and finds a
#Runge-Kutta methods with s steps and order d, using polynomial p
#in x and y
FindRKP:=proc(s,d,p,x,y) local B, a, b, drk, h, i, j, eq, var, sl:
    B:=GButcher(a, b, s):
    drk:=DiffRK(B, d, p, x, y, h);
    eq:={seq(coeff(drk, h, i)=0, i=0..d), add(b[i], i=1..s)=1}:
    var:={seq(b[i], i=1..s), seq(seq(a[i, j], j=1..i-1), i=1..s)}:
    sl:=solve(eq, var):
    B:=subs(sl, B):
    return B, sl:
end:

#FindRK(s,d): inputs positive integers s and d and finds all
#Runge-Kutta methods with s steps and order d.
FindRK:=proc(s,d) local x, y, c:
    return FindRKP(s, d, P(x,y,c,d+1), x, y)[1]:
end:

#FindRKR(s,d): inputs positive integers s and d and finds all
#Runge-Kutta methods with s steps and order d.
#Uses random polynomial for speed-up, may not be completely
#general/accurate
FindRKR:=proc(s,d) local x,y,M:
    M:=10:
    return FindRKP(s, d, RP(x,y,cap,d+1), x, y)[1]:
end:

#FindRK(s,d): inputs positive integers s and d and finds all
#Runge-Kutta methods with s steps and order d, smartly
#NOTE: this doesn't work, because for some reason you can't solve
#equations obtained from solve
FindRKSmart:=proc(s, d) local p, x, y, i, j, eq, var, sl, B:
    B:=GButcher(a, b, s):
    var:={seq(b[i], i=1..s), seq(seq(a[i, j], j=1..i-1), i=1..s)}:
    eq:={}:
    for i from 1 to d do
        p:=1+y^i:
        sl:=repair(FindRKP(s,d,p,x,y)[2]):
        eq:=solve(eq union sl, var):
        p:=1+x^i+y^i:
        sl:=repair(FindRKP(s,d,p,x,y)[2]):
        eq:=solve(eq union sl, var):
    od:
    B:=subs(eq, B):
    return B:
end: