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

Help:=proc():
    print(` CubicSplineCore(Lx,Lf,x) , CubicSplineF(Lx,Lf,x) , CubicSplineC(Lx,Lf,Lprimef,x) `):
end:

##PROBLEM 1##

#CubicSplineCore(Lx,Lf,x): Do the core work of both CubicSplineF
#and CubicSplineC, thereby avoiding the terrible practice of
#copying and pasting code
CubicSplineCore:=proc(Lx, Lf, x) local S, eq, var, a, i, j:
    if nops(Lx) <> 5 or nops(Lf) <> 5 then
        return FAIL:
    fi:
    #condition a
    S:=[seq(add(a[i,j]*x^j, j=0..3), i=0..3)]:
    var:={seq(seq(a[i,j], j=0..3), i=0..3)}:
    #condition b
    eq:={seq(subs(x=Lx[i], S[i])=Lf[i], i=1..4)}:
    eq:=eq union {subs(x=Lx[5], S[4])=Lf[5]}:
    #condition c
    eq:=eq union {seq(subs(x=Lx[i], S[i-1]) = 
        subs(x=Lx[i], S[i]), i=2..4)}:
    #condition d
    eq:=eq union {seq(subs(x=Lx[i], diff(S[i], x)) =
        subs(x=Lx[i], diff(S[i-1], x)), i=2..4)}:
    #condition e
    eq:=eq union {seq(subs(x=Lx[i], diff(S[i], x, x)) =
        subs(x=Lx[i], diff(S[i-1], x, x)), i=2..4)}:
    return S, eq, var:
end:

#CubicSplineF(Lx,Lf,x): inputs Lx and Lf, lists of numbers
#of length 5 with Lx increasing and x, a symbol for the variable
#outputs a list of cubic polynomials in x of length 4 that
#produces a cubic spline according to the specs at the bottom of
#p. 595, for the free and clamped boundary conditions,
#respectively. The ouptut should be a list of cubic polynomials
#of length 4,
#[S0,S1,S2,S3]
#according to the notation there.
CubicSplineF:=proc(Lx,Lf,x) local S, eq, var:
    S, eq, var:=CubicSplineCore(Lx, Lf, x):
    #condition fi
    eq:=eq union {subs(x=Lx[1], diff(S[1], x, x)) = 0,
        subs(x=Lx[4], diff(S[4], x, x)) = 0}:
    var:=solve(eq, var):
    return subs(var, S):
end:

#CubicSplineC(Lx,Lf,Lprimef,x) inputs Lx and Lf, lists of numbers
#of length 5 with Lx increasing, Lprimef, a list of length 2 that
#indicates the values of the derivate of the function at the
#endpoints x0 (alias Lx[1]), and x4 (alias Lx[5]),
#x, a symbol for the variable
#outputs a list of cubic polynomials in x of length 4 that
#produces a cubic spline according to the specs at the bottom of
#p. 595, for the free and clamped boundary conditions,
#respectively. The ouptut should be a list of cubic polynomials
#of length 4,
#[S0,S1,S2,S3]
#according to the notation there.
CubicSplineC:=proc(Lx,Lf,Lprimef,x) local S, eq, var:
    if nops(Lprimef) <> 2 then
        return FAIL:
    fi:
    S, eq, var:=CubicSplineCore(Lx, Lf, x):
    #condition fii
    eq:=eq union {subs(x=Lx[1], diff(S[1], x)) = Lprimef[1],
        subs(x=Lx[5], diff(S[4], x)) = Lprimef[2]}:
    var:=solve(eq, var):
    return subs(var, S):
end: