a = 0;
b = 5;

n = 100;

h = (b-a)/n;

x = a:h:b;

y = x.^2 + 1;

deriv = zeros(1, length(x)-1);

for i = 1:length(deriv)
    deriv(i) = (y(i+1) - y(i))/h;
end

plot(x(1:length(x)-1), deriv);