Initialization on non-DAE models

[ChrisRackauckas]

#2508 pointed out that the heuristic of implicit_dae || calculate_massmatrix(sys) !== I, i.e. "only do initialization on DAEs", while sensible for the DAE model, doesn't quite fit in the MTK context. Instead, what we need to do is check whether the varmap is correct and consistent for the chosen ODE. If it's not consistent for the selected ODE, then we still need to run an initialization step for the resulting ODE solve.

This needs a few changes in the ODE solver as though, since currently only DAE solvers will run the initialization step. Thus the test case:

using ModelingToolkit, OrdinaryDiffEq, Test
using ModelingToolkit: t_nounits as t, D_nounits as D

function System(;name)
    vars = @variables begin
        dx(t), [guess=0]
        ddx(t), [guess=0]
    end
    eqs = [
        D(dx) ~ ddx
        0 ~ ddx + dx + 1
    ]
    return ODESystem(eqs, t, vars, []; name)
end

@mtkbuild sys = System()
prob = ODEProblem(sys, [sys.dx => 1], (0,1)) # OK
prob = ODEProblem(sys, [sys.ddx => -2], (0,1), guesses = [sys.dx => 1])
sol = solve(prob, Rodas5P())
sol = solve(prob, Tsit5())

gives an erroneous success as it skips initialization, using prob.u0 == [0.0], the sentinel value since initialization is only run on DAEs. We will need to override that so that initialization is always run on any system that provides an initializeprob.

Fixes #2508

[ChrisRackauckas]

This also coincidentally adds steady state initializations as almost a free piece, and late binding initialization_eqs

[ChrisRackauckas]

Ran locally since CI is taken over by codecov bumps. Passes!