更好的时间循环结构只有一个计算时间步长的位置。

x_array = [x0]; y_array = [y0]; h = h_init;
  x = x0; y = y0;
  while x < x_end
    [y_new, err] = RK4_step_w_error(x,y,rhs,h);
    factor = abs(tol/err)^0.2;
    if factor >= 1
      y_array(end+1) = y = y_new;
      x_array(end+1) = x = x+h;
    end
    h = factor*h;
  end

对于代码中给定的数据

rhs = @(x,y) 20*cos(x);
  x0 = 0; y0 = 1; x_end = 6.5; tol = 1e-3; h_init = 1;

这给出了与精确解相对的结果

计算出的点正好位于精确解上，对于它们之间的线段，需要使用“密集输出”插值。或者作为第一个改进，仅包括来自半步计算的中间值。

function [ y_next, err] = RK4_step_w_error(x,y,rhs,h)
    y2 = RK4_step(x,y,rhs,h);
    y1 = RK4_step(x,y,rhs,h/2);
    y1 = RK4_step(x+h/2,y1,rhs,h/2);
    y_next = y1;
    err = (y2-y1)/15;
  end 
  
  function y_next = RK4_step(x,y,rhs,h)
    k1 = h*rhs(x,y);
    k2 = h*rhs(x+h/2,y+k1);
    k3 = h*rhs(x+h/2,y+k2);
    k4 = h*rhs(x+h,y+k3);
    y_next = y + (k1+2*k2+2*k3+k4)/6;
  end

修订版1

返回的误差是实际的步长误差。但是，步长控制所需的误差是单位步长误差或误差密度，即步长误差除以h

function [ y_next, err] = RK4_step_w_error(x,y,rhs,h)
    y2 = RK4_step(x,y,rhs,h);
    y1 = RK4_step(x,y,rhs,h/2);
    y1 = RK4_step(x+h/2,y1,rhs,h/2);
    y_next = y1;
    err = (y2-y1)/15/h;
  end

将示例改为在稳定平衡的两个分支之间振荡的简单双稳态模型

rhs = @(x,y) 3*y-y^3 + 3*cos(x);
  x0 = 0; y0 = 1; x_end = 13.5; tol = 5e-3; h_init = 5e-2;

给出了解、误差（针对ode45积分）和步长的曲线图

红色叉号表示被拒绝步骤的步长。

修订版2

函数值中的误差可用作5阶外推值的误差指导，使该方法成为外推模式下的5阶方法。由于它使用4阶误差预测5阶最佳步长，因此建议使用一个警告因子，代码在适当位置进行更改以

factor = 0.75*abs(tol/err)^0.2;

...

  function [ y_next, err] = RK4_step_w_error(x,y,rhs,h)
    y2 = RK4_step(x,y,rhs,h);
    y1 = RK4_step(x,y,rhs,h/2);
    y1 = RK4_step(x+h/2,y1,rhs,h/2);
    y_next = y1+(y1-y2)/15;
    err = (y1-y2)/15;
  end

在图中，步长适当增大，但误差显示更尖锐和更大的尖峰，该版本的方法显然不太稳定。