# Adaptation of existing code to multivariable system

I am trying to adapt the existing code from "Mathematica package for analysis and control of chaos in nonlinear system" by Gutierrez. This existing code calculates the Lyapunov Exponent by mapping method. I think it should work for more than 1D. It works for Logistic function (1D) but not for Henon function (2D).

 Orbit[map_,x0_,n_]:=NestList[map,x0,n];

LyapunovExp[map_, x0_, n_] :=
Plus @@ (Function[x, Evaluate[Log[Abs[D[map[x], x]]]]] /@
Drop[Orbit[map, x0, n + 500], 500])/n

Logistic[r_] := Function[x, r x (1 – x)];

Plot[lya = LyapunovExp[Logistic[r], 0.5, 500], {r, 0, 4},
PlotStyle -> Red]

Henon[r_] := Function[x, {1 - r x[[1]]^2 + x[[2]], 1/3 x[[1]]}]

Plot[LyapunovExp[Henon[r], {0.5, 0.5}, 150], {r, 0, 1.3},
PlotStyle -> Red]


My guess is that the problem starts with function D which is not a proper way to evaluate derivative of the multivariable function. Am I right?