**Short version**: My understanding of the use of TargetUnits in a plot statement is that if the function being traced has units, MMA will attempt to convert the output of that function to the unit form specified by TargetUnits. However, it does not seem to be doing that. What am I missing?

**Context**: I'm writing some general routines to analyze a particular type of differential equation system. These may be of arbitrary complexity, although with a well-defined relationship between units. I have not been able to get MMA to run NDSolve on systems that include units, so I have developed a way to convert all units to YES, remove the units, do NDSolve and then re-add the appropriate units to the resulting InterpolatingFunctions. When I look at the individual values that result from those interpolation functions, the results look good. But when I try to plot the results using TargetUnits, they do not.

**Simplified example**:

```
rateeq(t_) = {a'(t) == k1 b(t) + k2 a(t) b(t) + k3 b(t)^2,
b'(t) == -k1 b(t) - k2 a(t) b(t) - k3 b(t)^2};
initCond(t_) = {b(t) == Quantity(10^-9, "Molar"),
a(t) == Quantity(10^-5, "Molar")};
rateConstValues = {k1 -> Quantity(7000, 1/"Seconds"),
k3 -> Quantity(5.4 10^9, 1/("Molar" "Seconds")),
k2 -> Quantity(1.8 10^7, 1/("Molar" "Seconds"))};
quantitiesToSolveFor(t_) = {a(t), b(t)};
```

Here is my function to delete units after converting to SI:

```
UnitStrip(exp_) := Replace(exp, a_?QuantityQ :> QuantityMagnitude(UnitConvert(a)), All)
```

I've tried it and found it to work as expected. Here is my call to NDSolve:

```
result = NDSolve(UnitStrip(Join(rateeq(t),
initCond(Quantity(0, "Microseconds"))) /.
rateConstValues),
quantitiesToSolveFor(t),
UnitStrip({t, 0, Quantity(500, "Microseconds")}))((1))
```

This gives me replacement rules for a (t) and b (t) as interpolation functions. I can make a version of those same results but with the units put back as follows:

```
UnitizedResult = (# -> (Quantity(# /. result, "Moles"/"Meters"^3) /.
t -> (t/Quantity(1, "Seconds")))) & /@
(quantitiesToSolveFor(t))
```

Let's convert b (t) into a function and try a value:

```
f(t_) = b(t) /. UnitizedResult;
UnitConvert(f(Quantity(250, "Microseconds")), "Molar")
(* 1.65921 10^-10 M *)
```

And that result is correct. Here is the strange part:

```
Plot(Evaluate(f(t)),
{t, Quantity(0, "Microseconds"), Quantity(500, "Microseconds")},
TargetUnits -> "Molar")
```

Not only is the magnitude far apart (about 11 orders of magnitude), but the trend is wrong … it should be a decreasing function that approaches zero as an asymptote. I have verified that f (t) behaves as it should when using a manipulator:

```
Manipulate(UnitConvert(f(Quantity(t, "Microseconds")), "Molar"), {t, 0, 500})
```

Any ideas?

**Edited to add**: So, if I do not do the conversion on the t axis to a unified version, everything works fine. Does anyone know why my approach to transforming the t-axis is not working here? And does someone have a way to make the conversion successfully?