Evaluate expression in environment passed to function as parameter in R

I am trying to create a function for theoretical hessian matrix that I can then evaluate at different locations. First I tried setting expressions as values in a matrix or array, but although I could initially set an expression into a matrix I couldn't replace with the value calculated.

hessian_matrix <- function(gx, respect_to){
out_mat <- matrix(0, nrow=length(respect_to), ncol=length(respect_to))
  for(i in 1:length(respect_to)){
    for(j in 1:length(respect_to)){
      dthetad2x <- deriv(D(gx, respect_to[i]), respect_to[j], function.arg=TRUE)
      # also tried
      # dthetad2x <- as.expression(D(D(gx, respect_to[i])))
      out_mat[i,j] <- dthetad2x
    }
return(out_mat)
}

Because that didn't work, I decided to create an environment to house the indeces of the hessian matrix as object.

hessian_matrix <- function(gx, respect_to){
  out_env <- new.env()
  for(i in 1:length(respect_to)){
    for(j in 1:length(respect_to)){
      dthetad2x <- as.call(D(D(gx, respect_to[i]), respect_to[j]))
      assign(paste0(i,j), dthetad2x, out_env)
    }
  }
  return(out_env)
}

g <- expression(x^3-2*x*y-y^6)
h_g <- hessian_matrix(g, respect_to = c('x', 'y'))

This worked, and when I pass this in as a parameter to evaluate I can see the expression, but I can't evaluate it. I tried with call(), eval(), do.call(), get(), etc. and it didn't work. I also assigning the answer within the environment passed, making a new environment to return, or simply using variables.

fisher_observed <- function(h, at_val, params, sum=TRUE){
  out_env <- new.env()
  # add params to passed environment
  for(i in 1:length(at_val)){
    h[[names(at_val)[i]]] <- unname(at_val[i])
  }
  for(i in ls(h)){
    value <- do.call(i, envir=h, at_val)
    assign(i, value, out_env)
  }
  return(h)
}

fisher_observed(h_g, at_val=list(x=1,y=2))

According the code for do.call() this is how it should be used, but it isn't working when passed as a parameter in this way.

R already has the hessian matrix function. You do not have to write one. You could use deriv or deriv3 as shown below:

g <- expression(x^3 - 2 * x * y - y^6)
eval(deriv3(g, c('x','y')),list(x=1,y=2))

[1] -67
attr(,"gradient")
      x    y
[1,] -1 -194
attr(,"hessian")
, , x

     x  y
[1,] 6 -2

, , y

      x    y
[1,] -2 -480

If you want to use a function, you could do:

hessian <- function(expr,values){
  nms <- names(values)
  f <- eval(deriv3(g, nms),as.list(values))
  matrix(attr(f, 'hessian'), length(values), dimnames = list(nms,nms))
}

hessian(g, c(x=1,y=2))
   x    y
x  6   -2
y -2 -480

Although the function is not necessary as you would do double computation in case you wanted the gradient and hessian