function output = df_dx(f)

    eps = 1e-16;
    output = @(x) (f(x + eps) - f(x)) / eps;

end