# Manipulating a ParametricFunction to handle units

I have a pretty nice way to handle units on a result from NDSolve, which basically consists of converting all of the units to SI, stripping them, calling NDSolve, and then putting the units back. I'll show my example code below. I am trying, however, to update this routine to use ParametricNDSolve instead of NDSolve, and for the life of me I cannot figure out how to do something similar with the resulting ParametricFunctions. Any suggestions would be welcome.

First, my function for stripping units from any expression after converting them to SI units:

UnitStrip[exp_] := Replace[exp, a_?QuantityQ :> QuantityMagnitude[UnitConvert[a]], All]


Now my function for running NDSolve and adding the units back on:

UnitizedNDSolve[equations_, timeDependentVariables_,
initialConcentrations_, timeDependentVariableUnits_,
parameterValues_, parameterUnits_, {t_, t0_, tf_}] :=
Module[{initCond, result},
initCond = MapThread[(#1 /. t -> 0) == #2 &,
{timeDependentVariables, initialConcentrations}];
result = First@NDSolve[
UnitStrip[Join[equations, initCond] /. parameterValues],
timeDependentVariables, UnitStrip[{t, t0, tf}]];
QuantityUnit[UnitConvert[Quantity[1, #2]]]], #2] /.
{t -> (t/Quantity[1, "Seconds"])})) &,
{timeDependentVariables, timeDependentVariableUnits}]]


A test system:

eqns = {a'[t] == k b[t], b'[t] == -k b[t]};
tDepVars = {a[t], b[t]};
initConcentrations = {Quantity[1, "Molar"], Quantity[0.01, "Molar"]};
tDepVarUnits = QuantityUnit[initConcentrations];
paramVals = {k -> Quantity[6000, 1/"Seconds"]};
paramUnits = QuantityUnit[paramVals];


The function call and a plot of the results:

UnitizedResults = UnitizedNDSolve[eqns, tDepVars,
initConcentrations, tDepVarUnits,
paramVals, paramUnits,
{t, Quantity[0, "Seconds"], Quantity[500, "Microseconds"]}];
Plot[b[t] /. UnitizedResults /. t -> Quantity[time, "Microseconds"],
{time, 0, 500}, TargetUnits -> "Molar"] 