title: what is -nan(ind) and how do I resolve the error?
body: <p>Instead of getting numerical values, I get <code>-nan(ind)</code> when the number of iterations increases. I have applied the 4th order Runge Kutta method in the program, and the functions are given in <code>dvdt</code> and <code>drdt</code>.</p>
  <pre class="lang-cpp prettyprint-override"><code>#include &lt;iostream&gt;
  #include &lt;fstream&gt;
  #include &lt;cmath&gt;
  #include &lt;iomanip&gt;
  using namespace std;
  
  double dvdt(double t, double v);
  double drdt(double t, double r);
  
  double RKv(double t, double v, double h);
  double RKr(double t, double r, double h);
  
  double xo = 0;    //initial X position
  double yo = 0;    //initial Y position
  double x1;      //final X position
  double y11;      //final Y position
  double ti = 0;  //initial time = 0
  double tf = 30;      //time entered
  double dt = 0.01;   //time step (iteration time)
  double rud = 0;   //rudder angle
  double m = 53330750;   //mass of submarine (Kg)
  double Iz = 57.5;     //moment of inertia of submarine (kg.m^2)
  double xg = -0.012;   //x position of COG
  double yg = 0;     // position of COG
  double u = 10; // surge velocity
  double v, vo = 0; //sway velocity
  double r, ro = 0; // yaw rate
  double psio = 0; //initial ship heading
  double psi1;  //final ship heading
  double rudmax; // max rudder angle
  double au; //surge acceleration
  double av; // sway acceleration
  double ar; //angular acceleration
  double X; //surge force
  double Y; //sway force
  double N; //turning moment
  double h = 1;
  
  
  //matrices
  double M[2][2];         //mass matrix
  double Nu0[2][2];       //inertia matrix ig
  double b[2][1];         //force matrix
  double A[2][1];         //acceleration matrix
  double V[2][1];         //velocity matrix
  
  
  //maneuvering and hydrodynamic coefficents
  double Xu1 = -1.0467 * pow(10, -3);
  double Yv1 = -23.889 * pow(10, -3);
  double Yr1 = -1.0510 * pow(10, -3);
  double Nv1 = 1.1531 * pow(10, -3);
  double Nr1 = -0.50717 * pow(10, -3);
  double Xu = -2.8763 * pow(10, -3);
  double Yv = -38.948 * pow(10, -3);
  double Yr = 2.0031 * pow(10, -3);
  double Nv = -14.287 * pow(10, -3);
  double Nr = -4.2267 * pow(10, -3);
  double Yrud = 1.4308 * pow(10, -3);
  double Nrud = -0.71540 * pow(10, -3);
  
  double a1 = m - Yv1;
  double b1 = m * xg - Yr1;
  double a2 = m * xg - Yr1;
  double b2 = Iz - Nr1;
  double det = a1 * b2 - a2 * b1;
  
  double q1 = (-Yv * (Iz - Nr1) / det);
  double q2 = ((m - Xu1) * u - Yr) * (-m * xg - Yr1) / det;
  double q3 = ((Xu1 - Yv1) * u - Nv) * (-m * xg + Yr) / det;
  double q4 = ((m * xg - Yr1) * u - Nr)* (m - Yr1) / det;
  int main()
  {
      while (ti &lt;= tf)
      {
          vo = RKv(ti, vo, h);
          ro = RKr(ti, ro, h);
  
          v += vo;
          r += ro;
          
          cout &lt;&lt; &quot;v = &quot; &lt;&lt; v &lt;&lt; &quot;    &quot; &lt;&lt; &quot;r = &quot; &lt;&lt; r &lt;&lt; &quot;\n&quot;;
          ti = ti + 0.5;
          rud += 0.05;
      }
  }
  
  double RKv(double t, double v, double h)
  {
      double k1, k2, k3, k4, k5;
  
      
  
  
      for (int i = 0.1; i &lt;= h; i++)
      {
          
          k1 = h * dvdt(t, v);
          k2 = h * dvdt(t + 0.5 * h, v + 0.5 * k1);
          k3 = h * dvdt(t + 0.5 * h, v + 0.5 * k2);
          k4 = h * dvdt(t + h, v + k3);
  
          
          t = t + h;
          v = v + (1.0 / 6.0) * (k1 + 2 * k2 + 2 * k3 + k4);
  
          
      }
      return v;
  }
  
  double RKr(double t, double r, double h)
  {
      double k1, k2, k3, k4, k5;
  
      for (int i = 1; i &lt;= h; i++)
      {
          
          k1 = h * drdt(t, r);
          k2 = h * drdt(t + 0.5 * h, r + 0.5 * k1);
          k3 = h * drdt(t + 0.5 * h, r + 0.5 * k2);
          k4 = h * drdt(t + h, r + k3);
  
          
          t = t + h;
          r = r + (1.0 / 6.0) * (k1 + 2 * k2 + 2 * k3 + k4);
  
      }
      return r;
  }
  
  
  double dvdt(double t, double v)
  {
  //  return(-Yrud * rud - (-Yv * (Iz - Nr1) / det * v + ((m - Xu1) * u - Yr) * (-m * xg - Yr1) / det * r));
      return(-Yrud * rud - q1 * v + q2 * r);
      
  }
  
  double drdt(double t, double r)
  {
      //return(-Nrud * rud - (((Xu1 - Yv1) * u - Nv) * (-m * xg + Yr) / det * v + ((m * xg - Yr1) * u - Nr) * (m - Yr1) / det * r));
      return(-Nrud * rud - q3 * v + q4 * r);
      
  }
  </code></pre>
  <p>Here are the results of the program. I get numerical values for the first few iterations, after which I get <code>-nan(ind)</code>:</p>
  <p><a href="https://i.stack.imgur.com/wc8qN.png" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/wc8qN.png" alt="image" /></a></p>