compute_f = x -> [x[1]^2; x[2]^2 + x[3]; -x[3]^2 + x[4] + x[5]^3]

compute_Jf = x -> [
  2x[1] 0 0 0 0;
  0 2x[2] 1 0 0;
  0 0 -2x[3] 1 3*x[5]^2
]

function compute_Hf(x)
  Hf = zeros(3, 5, 5)

  Hf[1, 1, 1] = 2
  Hf[2, 2, 2] = 2
  Hf[3, 3, 3] = -2
  Hf[3, 5, 5] = 6 * x[5]

  return Hf
end

@testset "Tensors" begin
  m, n = 3, 5

  x = 2rand(n) .- 0.5

  @test TaylorTest.check(compute_f, compute_Jf, x)
  @test TaylorTest.check(compute_Jf, compute_Hf, x)

end