Я не понимаю, почему здесь происходит переполнение

Question

Извините за плохой английский.

Я старшекурсник по физике и плохо разбираюсь в программировании.

Я написал код для решения задачи нелинейного дифференциального уравнения по математике, и произошло переполнениеВот.Уравнения, которые я пытаюсь решить, следующие.х - позиция точечной массы, а М - номер точечной массы.бета просто константа.(Извините, я сейчас использую школьный компьютер, поэтому не могу скачать mathjax)

$ \ frac {d ^ 2} {dt ^ 2} x_ {i} = - x_ {i} - (x_{i} -x_ {i + 1}) - \ beta (x_ {i} ^ 3 + (x_ {i} - x {i + 1}) ^ 3) $, когда i = 1

$\ гидроразрыва {d ^ 2} {дт ^ 2} X_ {г} = - (X_ {я} -x_ {я-1}) - (X_ {я} -x_ {+ 1}) - \ бета ((x_ {i} - x_ {i-1}) ^ 3+ (x_ {i} -x_ {i + 1}) ^ 3) $, когда 2 <= i <= M-1 </p>

$ \гидроразрыва {d ^ 2} {дт ^ 2} X_ {М} = - (X_ {я} -x_ {я-1}) - X_ {я} -бета ((X_ {я} -x_ {я-1}) ^ 3- x_ {i} ^ 3) когда i = M

f[xvec_, i_, sizeM_] := 
  Which[i == 1, -2*xvec[[i]] + xvec[[i + 1]] - 
    beta*(xvec[[i]]^3 + (xvec[[i]] - xvec[[i + 1]])^3), i == sizeM, 
   xvec[[i - 1]] - 2*xvec[[i]] - 
    beta*((xvec[[i]] - xvec[[i - 1]])^3 - xvec[[i]]^3), 1 < i < sizeM,
    xvec[[i - 1]] - 2*xvec[[i]] + xvec[[i + 1]] - 
    beta*((xvec[[i]] - xvec[[i - 1]])^3 + (xvec[[i]] - 
          xvec[[i + 1]])^3)];
g[vvec_, i_] := vvec[[i]];

( f и g - уравнения, заданные )

sizeM = 16;(*the number of point mass*)
beta = 100.;

xvecInit = Table[Sin[N[Pi]*j/17.], {j, 1, sizeM}];

vvecInit = Table[0.0, {i, 1, sizeM}];
(*xvecInit and vvecInit are initial conditions*)

propT = 300.;(*time*)

storeT = 1.;
delT = 0.002;(*time interval*)

steps = 150000;(*time/time interval*)

Nstore = Floor[storeT/delT];

xvecOfT = Table[0.0, {j, 1, Floor[steps/Nstore] + 1}, {i, 1, sizeM}];
(*time evolution of the position*)
vvecOfT = Table[0.0, {j, 1, Floor[steps/Nstore] + 1}, {i, 1, sizeM}];
(*time evolution of the velocity*)

kx = Table[0.0, {j, 1, 4}, {i, 1, sizeM}];
kv = Table[0.0, {j, 1, 4}, {i, 1, sizeM}];
(*kx and kv are the vectors for using runge-kutta*)
xvecNow = xvecInit;
vvecNow = vvecInit;
xvecOfT[[1]] = xvecNow;
vvecOfT[[1]] = vvecNow;

    (* k1 *)
  Do[

    kx[[1, i]] = g[vvecNow, i];
    kv[[1, i]] = f[xvecNow, i, sizeM];

    , {i, 1, sizeM}];
  (* k2 *)
 Do[
    kx[[2, i]] = g[vvecNow + 0.5*delT*kv[[1]], i];
    kv[[2, i]] = f[xvecNow + 0.5*delT*kx[[1]], i, sizeM];
    , {i, 1, sizeM}];
  (* k3 *)
  Do[
   kx[[3, i]] = g[vvecNow + 0.5*delT*kv[[2]], i];
    kv[[3, i]] = f[xvecNow + 0.5*delT*kx[[2]], i, sizeM];
    , {i, 1, sizeM}];
  (* k4*)
  Do[
    kx[[4, i]] = g[vvecNow + delT*kv[[3]], i];
    kv[[4, i]] = f[xvecNow + delT*kx[[3]], i, sizeM];
    , {i, 1, sizeM}];

  xvecNow = 
    xvecNow + delT*(kx[[1]] + 2*kx[[2]] + 2*kx[[3]] + kx[[4]])/6.0;
  vvecNow = 
    vvecNow + delT*(kv[[1]] + 2*kv[[2]] + 2*kv[[3]] + kv[[4]])/6.0;

  If[Mod[k, Nstore] == 0,
    xvecOfT[[k/Nstore + 1]] = xvecNow;
    vvecOfT[[k/Nstore + 1]] = vvecNow;
    ];
  , {k, 1, steps}];

(*the problem happens in the last part. Overflow occurs and the computer cannot calculate*)

результатвыглядит следующим образом.

{{0.18375, 0.361242, 0.526432, 0.673696, 0.798017, 0.895163, 0.961826,
   0.995734, 0.995734, 0.961826, 0.895163, 0.798017, 0.673696, 
  0.526432, 0.361242, 0.18375}, {Overflow[], Overflow[], Overflow[], 
  Overflow[], Overflow[], Overflow[], Overflow[], Overflow[], 
  Overflow[], Overflow[], Overflow[], Overflow[], Overflow[], 
  Overflow[], Overflow[], Overflow[]}, {Overflow[], Overflow[], 
  Overflow[], Overflow[], Overflow[], Overflow[], Overflow[], 
  Overflow[], Overflow[], Overflow[], Overflow[], Overflow[], 
  Overflow[], Overflow[], Overflow[], Overflow[]}, {Overflow[], 
  Overflow[], Overflow[], Overflow[], Overflow[], Overflow[], 
  Overflow[], Overflow[], Overflow[], Overflow[], Overflow[], 
  Overflow[], Overflow[], Overflow[], Overflow[], 
  Overflow[]}, {Overflow[], Overflow[], Overflow[], Overflow[], 
  Overflow[], Overflow[], Overflow[], Overflow[], Overflow[],