title: Runge-Kutta 4th order method doesn't give the expected error
body: <p>I wrote a program to compare Euler and RK4 method in the simple problem of a ball thrown at a known angle and velocity with the only force of gravity. I calculated the absolute error and I graphed it with the time-step(h) to the x axis. The absolute error in Euler's turn out to be correct (proportional to h^2) but the absolute error of RK4 gives values from 0 to 10^-15. However I think the expected result is proportional to h^4. 
  Here is the code:</p>
  
  <pre><code>hr=0.1:0.05:1;h=0.1/4:0.05/4:1/4; %Chosen step for each method
  d1=zeros(1,19);d2=zeros(1,19); %difference between calculated and actual value for the two methods
  ne=0;nrk=0;
  for j=1:19%j counting the number of the step-sizes (19 different step-sizes)
  n=2500*hr(j);
  x=zeros(1,n);xrunge=zeros(1,n);
  y=zeros(1,n);yrunge=zeros(1,n);
  uy=zeros(1,n);uyrunge=zeros(1,n);
  t=zeros(1,n);tr=zeros(1,n);
  %Initializing
  uy(1)=20*sin(pi/3);ux=20*cos(pi/3);
  uyrunge(1)=20*sin(pi/3);
  a=-10;
  %Theoritical equations of motion for position and velocity to the y axis
  function [y]=f(t)
      y=20*sin(pi/3).*t-5*t.^2;
  end
  %Euler
  for i=2:n  
      if y(i-1)+h(j)*uy(i-1)&lt;0
          break;
      end
      if h(j)==0.1
          ne=ne+1;%number of calculations for uy with step-size 0.1
      end
      t(i)=t(i-1)+h(j);
      y(i)=y(i-1)+h(j)*uy(i-1);
      x(i)=x(i-1)+h(j)*ux;
      uy(i)=uy(i-1)+h(j)*a;  
      if i==2
          d1(j)=abs(y(i)-f(t(i))); %number of operations done to calculate y with step-size 0.1
      end
  end
  %RK4
  for i=2:n
      tr(i)=tr(i-1)+hr(j);
      xrunge(i)=xrunge(i-1)+hr(j)*ux;
      uyrunge(i)=uyrunge(i-1)+hr(j)*a;
      k1y=uyrunge(i-1);
      k2y=uyrunge(i-1)+hr(j)*a/2;
      k3y=uyrunge(i-1)+hr(j)*a/2;
      k4y=uyrunge(i-1)+hr(j)*a;
      if  yrunge(i-1)+hr(j)*(k1y+2*k2y+2*k3y+k4y)/6&lt;0
          break;
      end
      yrunge(i)=yrunge(i-1)+hr(j)*(k1y+2*k2y+2*k3y+k4y)/6;
      if i==2
          d2(j)=abs(yrunge(i)-f(tr(i))); %difference at each step-size for the first iteration
      end
  end
  end
  %function to fit |y-f(t)| according to least squares
  function [d]=f2(h)
  d=5*h.^2-5.621*10^(-7)*h+3.125*10^(-8);
  end
  figure(3)
  plot(h,d1,'r.',hr,d2,'k.',h,f2(h))
  xlabel('h')
  legend('Euler', 'RK4','Least squares fitting')
  ylabel( 'y-f(t)')
  title('Accuracy with step')
  disp(d2)
  </code></pre>
  
  <p>And the graph:
  <a href="https://i.stack.imgur.com/EBs02.png" rel="nofollow noreferrer">absolute error with time-step</a></p>