title: Comparison of odeint's runge_kutta4 with Matlab's ode45
body: <p>I would like to use <code>runge_kutta4</code> method in the <a href="http://headmyshoulder.github.io/odeint-v2/" rel="nofollow noreferrer">odeint C++ library</a>. I've solved the problem in Matlab. My following code in Matlab to solve <code>x'' = -x - g*x'</code>, with initial values <code>x1 = 1</code>, <code>x2 = 0</code>, is as follows</p>
  
  <p><code>main.m</code></p>
  
  <pre><code>clear all
  clc
  t = 0:0.1:10;
  x0 = [1; 0];
  [t, x] = ode45('ODESolver', t, x0);
  plot(t, x(:,1));
  title('Position');
  xlabel('time (sec)');
  ylabel('x(t)');
  </code></pre>
  
  <p><code>ODESolver.m</code></p>
  
  <pre><code>function dx = ODESolver(t, x)
  dx = zeros(2,1);
  g  = 0.15;
  dx(1) =  x(2);
  dx(2) = -x(1) - g*x(2);
  end
  </code></pre>
  
  <p>I've installed the odeint Library. My code for using <code>runge_kutta4</code> is as follows  </p>
  
  <pre><code>#include &lt;iostream&gt;
  
  #include &lt;boost/numeric/odeint.hpp&gt;
  
  using namespace std;
  using namespace boost::numeric::odeint;
  
  /* The type of container used to hold the state vector */
  typedef std::vector&lt; double &gt; state_type;
  
  const double gam = 0.15;
  
  /* The rhs of x' = f(x) */
  void lorenz( const state_type &amp;x , state_type &amp;dx ,  double t )
  {
      dx[0] =  x[1];
      dx[1] = -x[0] - gam*x[1];
  }
  
  
  int main(int argc, char **argv)
  {
      const double dt = 0.1;
      runge_kutta_dopri5&lt;state_type&gt; stepper;
      state_type x(2);
      x[0] = 1.0;
      x[1] = 0.0;
  
  
     double t = 0.0;
     cout &lt;&lt; x[0] &lt;&lt; endl;
     for ( size_t i(0); i &lt;= 100; ++i){
         stepper.do_step(lorenz, x , t, dt );
         t += dt;
         cout &lt;&lt; x[0] &lt;&lt; endl;
     }
  
  
      return 0;
  }
  </code></pre>
  
  <p>The result is in the following picture
  <img src="https://i.stack.imgur.com/T7tvP.png" alt="enter image description here"></p>
  
  <p>My question is why the result varies? Is there something wrong with my C++ code?</p>
  
  <p>These are the first values of both methods </p>
  
  <pre><code>Matlab    C++
  -----------------
  1.0000    0.9950
  0.9950    0.9803
  0.9803    0.9560
  0.9560    0.9226
  0.9226    0.8806
  0.8806    0.8304
  0.8304    0.7728
  0.7728    0.7084
  0.7083    0.6379
  </code></pre>
  
  <hr>
  
  <p>Update:</p>
  
  <p>The problem is that I forgot to include the initial value in my C++ code. Thanks for @horchler for noticing it. After including the proper values and using <code>runge_kutta_dopri5</code> as @horchler suggested, the result is </p>
  
  <pre><code>Matlab    C++
  -----------------
  1.0000    1.0000
  0.9950    0.9950
  0.9803    0.9803
  0.9560    0.9560
  0.9226    0.9226
  0.8806    0.8806
  0.8304    0.8304
  0.7728    0.7728
  0.7083    0.7084
  </code></pre>
  
  <p>I've updated the code to reflect these modifications. </p>