// *** 前進差分 ***
function df = diff_f1(x, dx, func)
   df = (func(x + dx) - func(x)) ./ dx
endfunction

// *** 中心差分 ***
function df = diff_f2(x, dx, func)
  df = (func(x + dx) - func(x - dx)) ./ (2 * dx)
endfunction

// *** 前進差分に対するRomberg1段 ***
function df = diff_f1r(x, dx, func)
  df = 2 * (diff_f1(x, dx, func) - 0.5 * diff_f1(x, 2 * dx, func))
endfunction

// *** 中心差分に対するRomberg1段 ***
function df = diff_f2r(x, dx, func)
  df = (4 / 3) * (diff_f2(x, dx, func) - 0.25 * diff_f2(x, 2 * dx, func))
endfunction