The answer to question two is ,"mostly yes" and it offers an almost immediate answer to your question:
funcD <- sapply(vars, function(v) D(func, v))
funcDD <- matrix(list(), 2,2)
for (i in 1:length(vars))
funcDD[,i] <- sapply(vars, function(v) D(funcD[[i]], v))
funcDD
#---------
[,1] [,2]
[1,] Expression Expression
[2,] Expression Expression
> funcDD[1,1]
[[1]]
-(sin(x + y) + cos(x - y))
The "mostly" qualification is that one needs to use "list" rather than "expression" as the object type that the matrix is holding. Expressions don't really qualify as atomic objects, and you could easily extract the value and use it as a call, which might even be more convenient than having it as an expression:
> is.expression(funcDD[1,1])
[1] FALSE
> funcDD[1,1][[1]]
-(sin(x + y) + cos(x - y))
> class(funcDD[1,1][[1]])
[1] "call"
Turns out what was wanted was the same structure, so this calls each matrix element with the same specific vector as the evaluation environment and returns them all as a matrix.:
matrix(sapply(funcDD, eval, env=list(x=0, y=pi)), length(vars))
#---------
[,1] [,2]
[1,] 1 -1
[2,] -1 1