onehot[dimensions_] := ArrayReshape[Flatten[IdentityMatrix[Times @@ dimensions]], Join[dimensions, dimensions]]
rotateTranspose[x_, k_] := TensorTranspose[x, RotateRight[Range[TensorRank[x]], k]];
NGrad[f_, x_, ϵ_] := 
rotateTranspose[Map[(((f[x + ϵ#]-f[x - ϵ#]))/(2ϵ)) &, onehot[Dimensions[x]], {TensorRank[x]}], TensorRank[x]];
NHessian[f_, x_, ϵ_] := Module[{
df=Function[X, Map[(((f[X + ϵ#]-f[X - ϵ#]))/(2ϵ)) &, onehot[Dimensions[x]], {TensorRank[x]}]]
},
rotateTranspose[Map[(((df[x + ϵ#]-df[x - ϵ#]))/(2ϵ))&, onehot[Dimensions[x]], {TensorRank[x]}], 2 TensorRank[x]]
]

(* scalar inputs/outputs *)
Π[x_, y_] := Sin[x+1] - Exp[x + y-1];

pΠpx[x_, y_]:=-(E^(-1+x+y))+Cos[1+x];
pΠpy[x_, y_]:=-(E^(-1+x+y));

p2Πpx2[x_, y_]:=-(E^(-1+x+y))-Sin[1+x];
p2Πpxpy[x_, y_]:=-(E^(-1+x+y));
p2Πpy2[x_, y_]:=-(E^(-1+x+y));

g[x_] := y0 /. First[NSolve[Π[x,y0]==0, y0, Reals]]

dgdx[x_] := Module[{y = g[x]},
-((pΠpy[x, y])^(-1)) pΠpx[x, y]
]

d2gdx2[x_] := Module[{y = g[x], dydx},
dydx = dgdx[x];
-(((pΠpy[x, y]))^(-1))((p2Πpxpy[x,y]+ p2Πpy2[x,y] dydx)dydx + (p2Πpx2[x,y]  + p2Πpxpy[x,y] dydx))
]

x0 = 0.5;
NGrad[g, x0, (10^(-6))] - dgdx[x0]

NHessian[g, x0, (10^(-4))] - NGrad[dgdx, x0, (10^(-4))]
NHessian[g, x0, (10^(-4))] - d2gdx2[x0]

Show[
ContourPlot[Π[x, y]==0, {x, -1, 2}, {y,-2, 2}],
Plot[g[x0] + dgdx[x0](x - x0)+(1/2)d2gdx2[x0](((x - x0))^2), {x, -1, 2}, PlotStyle->Red]
]

(* vector inputs/outputs *)
Π[{x1_, x2_}, {y1_, y2_}] := {x1 + y1 x2 - 0.5, Sin[x2 - x1 + (y2/y1)] + 0.4};

pΠpx[{x1_, x2_}, {y1_, y2_}]:={{1,y1},{-Cos[x1-x2-(y2/y1)],Cos[x1-x2-(y2/y1)]}};
pΠpy[{x1_, x2_}, {y1_, y2_}]:={{x2,0},{-((y2 Cos[x1-x2-(y2/y1)])/(y1^2)),((Cos[x1-x2-(y2/y1)])/y1)}};

p2Πpx2[{x1_, x2_}, {y1_, y2_}]:={{{0,0},{0,0}},{{Sin[x1-x2-(y2/y1)],-Sin[x1-x2-(y2/y1)]},{-Sin[x1-x2-(y2/y1)],Sin[x1-x2-(y2/y1)]}}};
p2Πpxpy[{x1_, x2_}, {y1_, y2_}]:={{{0,0},{1,0}},{{((y2 Sin[x1-x2-(y2/y1)])/(y1^2)),-((Sin[x1-x2-(y2/y1)])/y1)},{-((y2 Sin[x1-x2-(y2/y1)])/(y1^2)),((Sin[x1-x2-(y2/y1)])/y1)}}};
p2Πpy2[{x1_, x2_}, {y1_, y2_}]:={{{0,0},{0,0}},{{((2 y2 Cos[x1-x2-(y2/y1)])/(y1^3))+(((y2^2) Sin[x1-x2-(y2/y1)])/(y1^4)),-((Cos[x1-x2-(y2/y1)])/(y1^2))-((y2 Sin[x1-x2-(y2/y1)])/(y1^3))},{-((Cos[x1-x2-(y2/y1)])/(y1^2))-((y2 Sin[x1-x2-(y2/y1)])/(y1^3)),((Sin[x1-x2-(y2/y1)])/(y1^2))}}}

g[x_] := {y1, y2} /.FindRoot[Π[x, {y1, y2}], {{y1, 0.5}, {y2,0}}]

dgdx[x_] := Module[{y = g[x]},
-LinearSolve[pΠpy[x, y], pΠpx[x, y]]
]

d2gdx2[x_] := Module[{y = g[x], dydx},
dydx = dgdx[x];

output = ConstantArray[0, {2, 2, 2}];
Do[
output[[All, k, m]] = 
-LinearSolve[pΠpy[x, y],
p2Πpx2[x,y][[All,k,m]] + 
p2Πpxpy[x,y][[All,k]].dydx[[All, m]]+ 
p2Πpxpy[x,y][[All,m]].dydx[[All, k]]+
(p2Πpy2[x,y].dydx[[All, m]]).dydx[[All, k]]
]
,
{k, 1, 2}, {m, 1, 2}
];
output
]

{x1, x2} = {0.2, -0.4};
NGrad[g, {x1, x2}, (10^(-6))] - dgdx[{x1, x2}]
NGrad[dgdx, {x1, x2}, (10^(-6))]-d2gdx2[{x1, x2}]
NHessian[g, {x1, x2}, (10^(-5))]-d2gdx2[{x1, x2}]