Error in Lyapunov Exponent

I want to calculate Lyapunov exponents for by varying two parameters so that I can get density plot. I am using the code:

GramSchmidt[w_?MatrixQ] := 
 Module[{v = ConstantArray[0, Length[w]]}, 
  Table[v[[n]] = 
    w[[n]] - 
     Sum[(v[[i]] . w[[n]]/v[[i]] . v[[i]])*v[[i]], {i, n - 1}], {n, 
    Length[w]}];
  v]
LyapunovExponents[eqnsin_List, icsin : ({__Rule} | _Association), 
   nlein_Integer : 0, opts___?OptionQ] := 
  Module[{(*options*)tstep, maxsteps, ndsolveopts, logbase, showplot, 
    plotexponents, plotopts,(*other variables*)\[Delta], neq, nle, 
    vars, rhs, jac, eqns, unks, ics, cum, res, edat, state, newstate, 
    sol, W, norms},(*parse options*)
   tstep = Evaluate[TStep /. Flatten[{opts, Options[LyapunovExponents]}]];
   maxsteps = 
    Evaluate[MaxSteps /. Flatten[{opts, Options[LyapunovExponents]}]];
   ndsolveopts = 
    Evaluate[NDSolveOpts /. Flatten[{opts, Options[LyapunovExponents]}]];
   logbase = 
    Evaluate[LogBase /. Flatten[{opts, Options[LyapunovExponents]}]];
   showplot = 
    Evaluate[ShowPlot /. Flatten[{opts, Options[LyapunovExponents]}]];
   plotexponents = 
    Evaluate[
     PlotExponents /. Flatten[{opts, Options[LyapunovExponents]}]];
   plotopts = 
    Evaluate[PlotOpts /. Flatten[{opts, Options[LyapunovExponents]}]];
   neq = Length[eqnsin];
   If[nlein == 0, nle = neq, 
    nle = nlein];(*how many exponents*)(*extract vars and right hand s\
ides from eqnsin*)vars = eqnsin[[All, 1, 0, 1]];
   rhs = eqnsin[[All, 2]];
   (*jacobian matrix*)jac = D[rhs, {Replace[vars, {x_ -> x[t]}, 1]}];
   eqns = 
    Join[eqnsin, 
     Flatten[Table[\[Delta][i, j]'[
         t] == (jac . Table[\[Delta][i, j][t], {i, neq}])[[i]], {j, 
        nle}, {i, neq}]]];
   unks = Join[vars, Flatten[Table[\[Delta][i, j], {j, nle}, {i, neq}]]];
   ics = 
    Join[Table[var[0] == (var /. icsin), {var, vars}], 
     Flatten[Table[\[Delta][i, j][0] == 
        IdentityMatrix[neq][[i, j]], {j, nle}, {i, neq}]]];
   cum = Table[0, {nle}];
   state = 
    First@NDSolve`ProcessEquations[Flatten[Join[eqns, ics]], unks, t, 
      Evaluate[Sequence @@ ndsolveopts]];
   (*main loop*)
   edat = Table[newstate = First@NDSolve`Reinitialize[state, ics];
     NDSolve`Iterate[newstate, c tstep];
     sol = NDSolve`ProcessSolutions[newstate];
     W = 
      GramSchmidt[
       Evaluate[
        Table[\[Delta][i, j][c tstep], {j, nle}, {i, neq}] /. sol]];
     norms = Map[Norm, W];
     (*update running vector magnitudes*)
     cum = cum + Log[logbase, norms];
     ics = 
      Join[Table[var[c tstep] == (var[c tstep] /. sol), {var, vars}], 
       Flatten[Table[\[Delta][i, j][c tstep] == (W/norms)[[j, i]], {j,
           nle}, {i, neq}]]];
     cum/(c tstep), {c, maxsteps}];
   If[showplot, 
    Print[ListPlot[Transpose[edat][[plotexponents]], 
      Evaluate[Sequence @@ plotopts]]]];
   Return[cum/(maxsteps tstep)]];
Options[LyapunovExponents] = {NDSolveOpts -> {}, TStep -> 1, 
  MaxSteps -> 10^4(*10^4*), LogBase -> E, ShowPlot -> False, 
  PlotExponents -> 1(*{1,2,3}*), PlotOpts -> {}};

eqns1 = Table[a*x[t] - a*y[t] - y[t]*z[t], {a, 0, 2, 1}];
eqns2 = Table[-5* z[t] + d*x[t] + x[t]*y[t], {d, 0, 2, 1}];
eqns = {x'[t] == eqns1[[#]], 
   y'[t] == -1.67*y[t] + x[t]*z[t], {z'[t] == eqns2[[#]] & /@ 
     Range[1, Length[eqns2]]}} & /@ Range[1, Length[eqns1]]

ics = {x -> -1.0, y -> -1.0, z -> -1.0};
CA = Table[d, {d, 0, 2, 1}]
BA = Table[{a, CA}, {a, 0, 2, 1}]

AB = LyapunovExponents[eqns[[#]], ics, ShowPlot -> False] & /@ 
  Range[1, Length[eqns]]

LYa = {CA[[#]], BA[[#]], AB[[#, 1]], AB[[#, 2]], AB[[#, 3]]} & /@ 
  Range[1, Length[eqns]]

DensityPlot[LYa[[All, {1, 3}]][d, a], {d, 0, 2}, {a, 0, 2}, 
 PlotLegends -> Automatic]

I am not getting where I am doing mistakes. This code gives error.