Allow higher order derivatives in Boundary Conditions

[akashkgarg]

Currently the discretization assumes either Dirichlet boundary conditions or 1st order Neumann BCs. There are however PDEs that require higher order BCs; e.g., the beam equation with 3rd order BCs:

    @parameters x, t
    @variables u(..)
    Dt = Differential(t)
    Dtt = Differential(t)^2
    Dx = Differential(x)
    Dxx = Differential(x)^2
    Dx3 = Differential(x)^3
    Dx4 = Differential(x)^4

    g = -9.81
    EI = 1
    mu = 1
    L = 10.0
    dx = 0.4

    eq = Dtt(u(t,x)) ~ -mu*EI*Dx4(u(t,x)) + mu*g

    bcs = [u(0, x) ~ 0,
           u(t,0) ~ 0,
           Dx(u(t,0)) ~ 0,
           Dxx(u(t, L)) ~ 0,
           Dx3(u(t, L)) ~ 0]

    # Space and time domains
    domains = [t ∈ IntervalDomain(0.0,1.0),
               x ∈ IntervalDomain(0.0,L)]

    pdesys = PDESystem(eq,bcs,domains,[t,x],[u(t,x)])
    discretization = MOLFiniteDifference([x=>dx],t, centered_order=4)
    prob = discretize(pdesys,discretization)

A related issue is also reported here with second order BCs:
https://github.com/SciML/DiffEqOperators.jl/issues/344