Array Function Generation for MOL

[ChrisRackauckas]

We discussed that we can get more optimized codegen here if we target the array form:

D(u[1]) ~ u[1]
D(u[2:end-1]) ~ u[1:n-2] .- 2u[2:end-1] .+ u[3:end] .+ ifelse(u[2:end-1]>0, a, b)*u[2:end-1]
D(u[end]) ~ u[n] - u[n-1]

Which would then hit the broadcast CSE in the MTK codegen and thus give optimized looping code. We can then define behind the scenes some observed equations u(x, y, ⋅) ~ u[1:n,1:n] which would then make sol[u(x, y, ⋅)] -> sol[u[1:n,1:n]] and the specialized PDE interfaces shouldn't be required for the solution handling.