Решить систему дифференциальных уравнений

Я пытаюсь решить систему дифф. уравнений в Wolfram Mathematica, используя метод Рунге-Кутты. Необходимо найти решения на промежутке [0, 100], но в какой-то момент программа перестает находить корни системы. Что я делаю не так?

Также нужно найти погрешности, и я сравнивал полученные решения с решением NDSolve, и оно дает слишком большие значения.

tet = 10;
r = 28;
b = 0.1;
e = (r - 1)^(-1/2);
a = e*b*tet^(-1/2);
lmb = e*(b*tet + 1)^(-1/2);
Print["e=", e]
Print["a=", a]
Print["lmb=", lmb]
xk = 2;
yk = 1;
zk = 1;
h = 0.1;
fx [t_, x_, y_, z_] = y;
fy [t_, x_, y_, z_] = x - lmb*y - x*z;
fz [t_, x_, y_, z_] = -a*z + x^2;
solution = 
  NDSolve[{x'[t] == y[t], y'[t] == x[t] - lmb*y[t] - x[t]*z[t], 
    z'[t] == -a*z[t] + x[t]^2, x[0] == 2, y[0] == 1, z[0] == 1}, {x, 
    y, z}, {t, 100}];
Do[{Print["Шаг h=", h]; 
  res = NSolve[{
     k1x == 
      h*fx[t + (1/2 - Sqrt[3]/6) h, 
        xk + 1/4 k1x + (1/4 - Sqrt[3]/6) k2x, 
        yk + 1/4 k1y + (1/4 - Sqrt[3]/6) k2y, 
        zk + 1/4 k1z + (1/4 - Sqrt[3]/6) k2z],
        k2x == 
      h*fx[t + (1/2 + Sqrt[3]/6) h, 
        xk + (1/4 + Sqrt[3]/6) k1x + 1/4 k2x, 
        yk + (1/4 + Sqrt[3]/6) k1y + 1/4 k2y, 
        zk + (1/4 + Sqrt[3]/6) k1z + 1/4 k2z],
     k1y == 
      h*fy[t + (1/2 - Sqrt[3]/6) h,  
        xk + 1/4 k1x + (1/4 - Sqrt[3]/6) k2x, 
        yk + 1/4 k1y + (1/4 - Sqrt[3]/6) k2y, 
        zk + 1/4 k1z + (1/4 - Sqrt[3]/6) k2z],
        k2y == 
      h*fy[t + (1/2 + Sqrt[3]/6) h, 
        xk + (1/4 + Sqrt[3]/6) k1x + 1/4 k2x, 
        yk + (1/4 + Sqrt[3]/6) k1y + 1/4 k2y, 
        zk + (1/4 + Sqrt[3]/6) k1z + 1/4 k2z],
     k1z == 
      h*fz[t + (1/2 - Sqrt[3]/6) h,  
        xk + 1/4 k1x + (1/4 - Sqrt[3]/6) k2x, 
        yk + 1/4 k1y + (1/4 - Sqrt[3]/6) k2y, 
        zk + 1/4 k1z + (1/4 - Sqrt[3]/6) k2z],
        k2z == 
      h*fz[t + (1/2 + Sqrt[3]/6) h, 
        xk + (1/4 + Sqrt[3]/6) k1x + 1/4 k2x, 
        yk + (1/4 + Sqrt[3]/6) k1y + 1/4 k2y, 
        zk + (1/4 + Sqrt[3]/6) k1z + 1/4 k2z]}, 
    {k1x, k2x, k1y, k2y, k1z, k2z}, Reals]  ;
  xk1 = xk + h/2 (res[[1, 1, 2]] + res[[1, 2, 2]]);
  yk1 = yk + h/2 (res[[1, 3, 2]] + res[[1, 4, 2]]);
  zk1 = zk + h/2 (res[[1, 5, 2]] + res[[1, 6, 2]]);
  Print["xk: ", xk1];
  Print["yk: ", yk1];
  Print["zk: ", zk1];
  xk = xk1;
  yk = xk1;
  zk = xk1;
  solution = 
   NDSolve[{x'[t] == y[t], y'[t] == x[t] - lmb*y[t] - x[t]*z[t], 
     z'[t] == -a*z[t] + x[t]^2, x[0] == 2, y[0] == 1, z[0] == 1}, {x, 
     y, z}, {t, 100}];
  realx = x[h] /. solution;
  realy = y[h] /. solution;
  realz = z[h] /. solution;
  Print["real x=", realx];
  Print["real y=", realy];
  Print["real z=", realz];
  Print["Погрешность для x=", Abs[realx - xk1]]
   Print["Погрешность для y=", Abs[realy - yk1]]
   Print["Погрешность для z=", Abs[realz - zk1]]
  }, {h, 0, 1, 0.1}]

Я пробовал использовать Solve, NSolve, FindRoot но ничего не выдает корректных результатов.