struct MOLFiniteDifference{G} <: DiffEqBase.AbstractDiscretization

# Constructors. If no order is specified, both upwind and centered differences will be 2nd order
function MOLFiniteDifference(dxs, time=nothing; approx_order = 2, upwind_order = 1, grid_align=CenterAlignedGrid())
    if approx_order % 2 != 0
        @warn "Discretization approx_order must be even, rounding up to $(approx_order+1)"
    @assert approx_order >= 1 "approx_order must be at least 1"
    @assert upwind_order >= 1 "upwind_order must be at least 1"
    return MOLFiniteDifference{typeof(grid_align)}(dxs, time, approx_order, upwind_order, grid_align)
eq = [your system of equations, see examples for possibilities]
bcs = [your boundary conditions, see examples for possibilities]

domain = [your domain, a vector of Intervals i.e. x ∈ Interval(x_min, x_max)]

@named pdesys = PDESystem(eq, bcs, domains, [t, x, y], [u(t, x, y)])

discretization = MOLFiniteDifference(dxs, 
                                      <your choice of continuous variable, usually time>; 
                                      advection_scheme = <UpwindScheme() or WENOScheme()>, 
                                      approx_order = <Order of derivative approximation, starting from 2> 
                                      grid_align = <your grid type choice>)
prob = discretize(pdesys, discretization)

Where dxs is a vector of pairs of parameters to the grid step in this dimension, i.e. [x=>0.2, y=>0.1]. For a non uniform rectilinear grid, replace any or all of the step sizes with the grid you'd like to use with that variable, must be an AbstractVector but not a StepRangeLen.

Note that the second argument to MOLFiniteDifference is optional, all parameters can be discretized if all required boundary conditions are specified.

Currently implemented advection schemes are UpwindScheme() and WENOScheme(), defaults to upwind. See advection schemes for more information.

Currently supported grid types: center_align and edge_align. Edge align will give better accuracy with Neumann boundary conditions.

center_align: naive grid, starting from lower boundary, ending on upper boundary with step of dx

edge_align: offset grid, set halfway between the points that would be generated with center_align, with extra points at either end that are above and below the supremum and infimum by dx/2. This improves accuracy for Neumann BCs.

Any unrecognized keyword arguments will be passed to the ODEProblem constructor, see its documentation for available options.