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

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

Help:=proc():
    print(` Zg(f, x, y, y0, N) , PolAppx(f,x,y,y0,x1,N,eps) `):
end:

##PROBLEM 1##
 
#Zg(f, x, y, y0, N): the first N+1 Maclaurin polynomials
#of solutions of y'(x)=f(x,y(x)), y(0)=y0
#Now, f can be any analytic function of x and y
Zg:=proc(f, x, y, y0, N) local L, n, cur:
    if diff(f,y)=0 then
        print(`the function f only depends on x, so the solution is the indefinite integral`):
        return y0+int(f,x=0..x):
    fi:
    
    L:=[y0]:
    cur:=y0:
    for n from 1 to N do
        cur:=Z1(mtaylor(f, [x, y], n+1), x, y, cur):
        L:=[op(L), cur]:
    od:
    return L:
end:

#Zg(cos(x*y)+exp(x),x,y,1,10)[11] = 1+2*x+1/2*x^2+1/6*x^3-11/24*x^4-59/120*x^5-139/720*x^6+1/5040*x^7+5353/13440*x^8+129473/362880*x^9+529633/3628800*x^10
#Zg(exp(sin(x)+y)+1+x+y,x,y,1,10)[11] = 1+4*x+15/2*x^2+205/12*x^3+3559/72*x^4+1332071/8640*x^5+1539270473/3110400*x^6+25361870561021/15676416000*x^7+169730206005400697/31603654656000*x^8+207237163591913418208537/11468334201569280000*x^9+256007310396822520065313781353/4161629115065460326400000*x^10
#Zg(cos(x*y)+exp(x)+1,x,y,1,10)[11] = 1+3*x+1/2*x^2+1/6*x^3-17/24*x^4-119/120*x^5-199/720*x^6+61/5040*x^7+13963/13440*x^8+379457/362880*x^9+1313569/3628800*x^10

#NOTE: taylor (1 variable) didn't work with the second one, so I used mtaylor for this problem

##PROBLEM 2##

#PolAppx(f,x,y,y0,x1,N,eps):  that inputs
#(i) an expression f in x and y that is analytic in x and y near the origin (0,0)
#(ii) variable names x,y
#(iii) a number y0
#(iv) a number x1
#(v) an integer N
#(vi) a small positive number eps
#and first checks that abs(x1) is less than the
#ERC(Zg(f,x,y,y0,N)/2, x)[1] (to be safe!), and if does not,
#returns FAIL. Then finds the MacLaurin polynomial of least
#degree (degree d) such that
#abs(subs(x=x1, Zg(f,x,y,y0,N)[d+1] - Zg(f,x,y,y0,N)[d])) < eps
#followed by its value at x=x1.
#If N was too small, returns FAIL
#NOTE: if you want eps/100, divide your own eps by 100. I'm not
#going to do that for you.
PolAppx:=proc(f, x, y, y0, x1, N, eps) local ZGL, d:
    ZGL:=Zg(f, x, y, y0, N):
    if abs(x1) >= ERC(ZGL[-1]/2, x)[1] then
        #No convergence here, sorry
        print(`No convergence here`):
        return FAIL:
    fi:
    for d from 1 to nops(ZGL) - 1 do
        if abs(subs(x=x1, ZGL[d+1] - ZGL[d])) < eps then
            return ZGL[d], subs(x=x1, ZGL[d]):
        fi:
    od:
    #Didn't converge fast enough
    #I'm not sorry this time
    print(`Didn't converge fast enough`):
    return FAIL:
end:

#NOTE: I changed all the 10^(-10)'s to 10^(-12)'s, since I'm
#not dividing eps by 100
#PolAppx(1+x^2+y^3,x,y,1,0.3,50,10^(-12)) returns 1+x, 1.3
#PolAppx(cos(x)*exp(y)+1+x+y,x,y,1,0.3,50,10^(-12)) runs for a really long time and probably returns FAIL (all the way through N=30 did).  The fail is the "Didn't converge fast enough" one.  This makes sense, because the radius of convergence is about 0.31.

##PROBLEM 3##

#See hw12.pdf