ref and x0 are scalar, and f produces a scalar when its input is a scalar.
Inside the for loop, h is a scalar, hence derivative is scalar as well. Thus, error becomes a scalar, which means that the line
rate = [error(1); error(2:end)./error(1:end-1)];
just evaluates to
rate = error;
and causes the behaviour you see.
There are too many things going wrong here to give you the fix with minimum changes, but this would probably give you what you need:
errors = [];
derivatives = [];
for h = hs
derivative = (f(x0+h) - f(x0-h))/(2*h);
derivatives = [derivatives; derivative];
err = abs(ref - derivative); % try not to shadow the default error()
errors = [errors; err];
end
rate = [errors(1); errors(2:end)./errors(1:end-1)];
disp([hs', derivatives, errors, rate])
