%  This is the cycle of the Lorenz_96 model

   %Runge-Kutta coefficients
   k1=timestep*deriv(forcing,truth(indices(1,:)),truth(indices(2,:)),truth(indices(3,:)),truth(index));
   k2=timestep*deriv(forcing,truth(indices(1,:))+k1(indices(1,:))/2.,truth(indices(2,:))+k1(indices(2,:))/2.,truth(indices(3,:))+k1(indices(3,:))/2.,truth(index)+k1(index)/2.);
   k3=timestep*deriv(forcing,truth(indices(1,:))+k2(indices(1,:))/2.,truth(indices(2,:))+k2(indices(2,:))/2.,truth(indices(3,:))+k2(indices(3,:))/2.,truth(index)+k2(index)/2.);
   k4=timestep*deriv(forcing,truth(indices(1,:))+k3(indices(1,:)),truth(indices(2,:))+k3(indices(2,:)),truth(indices(3,:))+k3(indices(3,:)),truth(index)+k3(index));

   truth(:)=truth(:)+(k1+2.*k2+2.*k3+k4)/6.;

   for j=1:ensemblesize
     k1=timestep*deriv(forcing,ensemble(indices(1,:),j),ensemble(indices(2,:),j),ensemble(indices(3,:),j),ensemble(index,j));
     k2=timestep*deriv(forcing,ensemble(indices(1,:),j)+k1(indices(1,:))/2.,ensemble(indices(2,:),j)+k1(indices(2,:))/2.,ensemble(indices(3,:),j)+k1(indices(3,:))/2.,ensemble(index,j)+k1(index)/2.);
     k3=timestep*deriv(forcing,ensemble(indices(1,:),j)+k2(indices(1,:))/2.,ensemble(indices(2,:),j)+k2(indices(2,:))/2.,ensemble(indices(3,:),j)+k2(indices(3,:))/2.,ensemble(index,j)+k2(index)/2.);
     k4=timestep*deriv(forcing,ensemble(indices(1,:),j)+k3(indices(1,:)),ensemble(indices(2,:),j)+k3(indices(2,:)),ensemble(indices(3,:),j)+k3(indices(3,:)),ensemble(index,j)+k3(index));

     ensemble(:,j)=ensemble(:,j)+(k1+2.*k2+2.*k3+k4)/6.;
   end