AutoDiff

AutoDefineDerivatives[FnNameString_String,PureFunExpression_,VectorLength_Integer,MaxDer_Integer]:= Module[
{fOut=ToExpression[FnNameString],fTemp, Id=IdentityMatrix[VectorLength]},
fTemp[0]=PureFunExpression[#]& ;
Print["Redefining variable ", FnNameString];
Off[Part::"partd"];
 Do[
fTemp[order+1]=Map[Derivative[#][fTemp[order]]&,Id];
fTemp[order+1]=Function[x,Evaluate[Table[fTemp[order+1][[i]][x],{i,VectorLength}]]];,
{order,0,MaxDer-1}]; 
fOut["Count"]=ConstantArray[0,MaxDer+1];
Do[
fOut["Data",order]=fTemp[order],
{order,0,MaxDer}
];
fOut[order_]:=Module[{n=Length[fOut["Count"]]},
fOut["Count"]+=UnitVector[n,order+1];
(* Sow[#,{fOut,order}]&; FIX *)
fOut["Data",order]
];
On[Part::"partd"];
Print["Defined self counting fn and derivatives for ", FnNameString];
Print[FnNameString,"[\"Count\"] gives current count"];
Print[FnNameString,"[i] returns ith derivative and increments counters. I should make it sow pts"];
Print[FnNameString,"[\"Data\",i] gives ith derivative without the counter"];
]

Clear[MyFun,f,df]
f[{x1_,x2_,x3_}]:= Sin[x2 x3 + x1];
g=f;
dg[{x1_,x2_,x3_}]={Derivative[{1,0,0}][g][{x1,x2,x3}],Derivative[{0,1,0}][g][{x1,x2,x3}],Derivative[{0,0,1}][f][{x1,x2,x3}]};
ddg[{x1_,x2_,x3_}]={Derivative[{1,0,0}][dg][{x1,x2,x3}],Derivative[{0,1,0}][dg][{x1,x2,x3}],Derivative[{0,0,1}][dg][{x1,x2,x3}]};
AutoDefineDerivatives["f",f,3,2]

xx={1.1,2.1,1.3};
f["Count"]
{
{g[xx]-f[0][xx]},
{dg[xx]-f[1][xx]},
{ddg[xx]-f[2][xx]},
{ddg[xx]-f[2][xx]}
}
f["Count"]