Fix functions with an obligatory unitful domain

src/UnitfulRecipes.jl

@@ -76,11 +76,30 @@ end
     x, y, f.(x',y)
 end
 @recipe function f(f::Function, u::Units)
-    recipedata = RecipesBase.apply_recipe(plotattributes,f)
-    (f, xmin, xmax) = recipedata[1].args
-    f, xmin*u, xmax*u
+    uf = UnitFunction(f, [u])
+    recipedata = RecipesBase.apply_recipe(plotattributes, uf)
+    (_, xmin, xmax) = recipedata[1].args
+    return f, xmin*u, xmax*u
 end
 
+"""
+```julia
+UnitFunction
+```
+A function, bundled with the assumed units of each of its inputs.
+
+```julia
+f(x, y) = x^2 + y
+uf = UnitFunction(f, u"m", u"m^2")
+uf(3, 2) == f(3u"m", 2u"m"^2) == 7u"m^2"
+```
+"""
+struct UnitFunction <: Function
+    f::Function
+    u::Vector{Units}
+end
+(f::UnitFunction)(args...) = f.f((args .* f.u)...)
+
 #===============
 Attribute fixing
 ===============#

test/runtests.jl

@@ -122,6 +122,8 @@ end
         pl = plot(f, m)
         @test xguide(pl) == string(m)
         @test yguide(pl) == string(m^2)
+        f(x) = exp(x/(3m))
+        @test plot(f, u"m") isa Plots.Plot
     end
 end