This works:
julia> prob = ODEProblem(ode3a, [300.0u"K"], tspan, params)
ODEProblem with uType Vector{Quantity{Float64, π―, Unitful.FreeUnits{(K,), π―, nothing}}} and tType Quantity{Float64, π, Unitful.FreeUnits{(s,), π, nothing}}. In-place: true
timespan: (0.0 s, 6.0 s)
u0: 1-element Vector{Quantity{Float64, π―, Unitful.FreeUnits{(K,), π―, nothing}}}:
300.0 K
julia> solve(prob, Tsit5())
retcode: Success
Interpolation: specialized 4th order "free" interpolation
t: 5-element Vector{Quantity{Float64, π, Unitful.FreeUnits{(s,), π, nothing}}}:
0.0 s
...
Note the 300.0 and the Tsit5. The former is necessary to get the IC in the right format. The latter is using an explicit solver (specifying none seems to select Rodas5). This makes it much easier as there is no autodiff happening, which again adds constraints on the types.
It can still be done, but gets a bit tricker. Potentially easier, if you need an implicit solver, is to set autodiff=false, although this still throws an error:
julia> solve(prob, Rodas5(autodiff=false))
ERROR: DimensionError: K and false are not dimensionally compatible.
Stacktrace:
[1] convert(#unused#::Type{Quantity{Float64, π―, Unitful.FreeUnits{(K,), π―, nothing}}}, x::Bool)
not sure whatβs the problem there.