diff --git a/Project.toml b/Project.toml
index a505891..7c4e462 100644
--- a/Project.toml
+++ b/Project.toml
@@ -16,6 +16,7 @@ Test = "1.11.0"
 
 [extras]
 SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
+TensorOperations = "6aa20fa7-93e2-5fca-9bc0-fbd0db3c71a2"
 
 [targets]
-test = ["SpecialFunctions"]
+test = ["SpecialFunctions", "TensorOperations"]
diff --git a/test/tensors.jl b/test/tensors.jl
index ae4f8b7..47e4406 100644
--- a/test/tensors.jl
+++ b/test/tensors.jl
@@ -26,3 +26,74 @@ end
   @test TaylorTest.check(compute_Jf, compute_Hf, x)
 
 end
+
+TP = TO.tensorproduct
+
+function check_symmetry(a)
+  R = true
+  for i in axes(a, 1)
+    r = a[i,:,:]' == a[i,:,:]
+    R &= r
+    if !r
+      n = maximum(abs, a[i,:,:]' - a[i,:,:])
+      @warn "Tensor not symmetric at i=$i: $n"
+    end
+  end
+  return R
+end
+
+function prod(A, B)
+  return TP([1, 2, 3], A, [1, 2], B, [3])
+end
+
+function revprod(A, B)
+  return TP([1, 3, 2], B, [3], A, [1, 2])
+end
+
+function symprod(A, B)
+  return prod(A, B) + revprod(A, B)
+end
+
+@testset "Product of scalar and vector functions" begin
+  v = x -> [exp(2 * x[1]), cos(x[1] * x[2]), 1.0]
+  Jv = x -> [
+    2*exp(2 * x[1]) 0;
+    -x[2]*sin(x[1] * x[2]) -x[1]*sin(x[1] * x[2]);
+    0 0
+  ]
+  Hv = x -> begin
+    h = zeros(3, 2, 2)
+    h[1, :, :] = [
+      4*exp(2 * x[1]) 0;
+      0 0
+    ]
+
+    h[2, :, :] = [
+      -x[2]^2*cos(x[1] * x[2]) -sin(x[1] * x[2])-x[1]*x[2]*cos(x[1] * x[2]);
+      -sin(x[1] * x[2])-x[1]*x[2]*cos(x[1] * x[2]) -x[1]^2*cos(x[1] * x[2])
+    ]
+    h
+  end
+
+  φ = x -> x[1] * x[2]^2
+  ∇φ = x -> [x[2]^2; 2x[1]*x[2]]
+  Hφ = x -> [
+    0 2x[2];
+    2x[2] 2*x[1]
+  ]
+
+  f = x -> v(x) * φ(x)
+  Jf = x -> φ(x) * Jv(x) + v(x)*∇φ(x)'
+  Hf = x -> φ(x) * Hv(x) + symprod(Jv(x), ∇φ(x)) + TP(v(x), [1], Hφ(x), [2, 3])
+
+  x = 2rand(2) .- 0.5
+
+  @test TaylorTest.check(v, Jv, x)
+  @test TaylorTest.check(Jv, Hv, x)
+
+  @test TaylorTest.check(φ, ∇φ, x)
+  @test TaylorTest.check(∇φ, Hφ, x)
+
+  @test TaylorTest.check(f, Jf, x)
+  @test TaylorTest.check(Jf, Hf, x)
+end