Error $L_{2}$ convergence in Finite Element for Poisson Equation

Question

I have written a Matlab code to solve the equation $-u'' = f$ with conditions $u(0) = u'(1) = 0$ on the domain $x \in [0,1]$. I tested the code with $f(x) = -2, \forall x\in [0,1]$. I check the plot with one of element $n=8$, and the result $u = u^{h}$ agrees:

However, I am testing the convergence error $L_{2}$ for the problem $n=8,16,32,64$. The error plot I get is not giving me the result of $O(\Delta x^{2})$ where $\Delta x$ is the fixed step size. I am trying to find if there is a bug. Why is a big change in my error between $n=16$ and $n=32$? It is probably something small in my code, but for some reason I can't figure it out. Thanks in advance!!!! Furthermore, I use a uniform grid $[0,1]$ so $\Delta x (h)$ is a fixed step size.

This is the code that I solve with Matlab and a driver to test the convergence:

% script to generate the uniform mesh for Finite element: 
% dom = [a b] where a < b is the domain of our problem 
function [xgrid, h] = generate1d_uniform(n) 
h=1/(n-1); 
xgrid=linspace(0,1,n);
end 

% implement the solver for Finite Element for Poisson equation
function uh = fem1d(n,f)

% generate the uniform mesh: 
[x,h] = generate1d_uniform(n); 
% load vector 
b = loadVector1D(x,h,f); 
A = StiffMat1D(x,h); 

% since u0 is 0 
Ae=A(2:end,2:end); 
fe=b(2:end); 
uhe=Ae\fe; 

% adjust the boundary condition u(0)=0
uh=[0;uhe]; 
end 

% function to assembly load vector b: 
% f is the source function; x is from the domain generate by 1d mesh
function b = loadVector1D(x,h,f)
n = length(x)-1; 
b = zeros(n+1,1); 

for i = 1:n
   b(i) = b(i) + f(x(i))*(h/2); 
   b(i+1) = b(i+1) + f(x(i+1))*(h/2); 
end

end 

% function to load the stiffness matrix for A: 
% the local stiffness matrix is of the form Ak = 1/h [1 -1; -1 1] 
% adjust the A(1,1) = 1 and A(n+1,n+1) = 1 for boundary conditions in this
% case u'(1) = 0 and u(0) = 0 
function A = StiffMat1D(x,h)
n = length(x)-1; 
Ak = spdiags(ones(n+1,1)*[-1 2 -1],-1:1,n+1,n+1); 
% Adjust the boundary condition u'(1) = 0
Ak(1,1) = 1; 
Ak(end,end) = 1;   
A = (1/h)*Ak; 
end

%%% DRIVER for convergence error %%%
function driver1
clear all; 
clc;
format short; 
%number of elements set: 
nset=[8,16,32,64]; 

% source function in part (a) 
f = @(x) -2.*(0<=x & x<=1);

% exact solution: 
u = @(x) x.^2 - 2*x; 
%define the error L2 set to store the values: 
errorL2=zeros(size(nset)); 
hset=zeros(size(nset));
% iteration: 
niter=length(nset); 

    for i=1:niter
        n=nset(i); 
        disp(n)
        errorL2(i)=compute_error(u,n,f);
    end

% write the table for error L2: 
T1=table(nset', errorL2'); 
T1.Properties.VariableNames ={'n elements', 'Error L2'}; 
writetable(T1,'TableError.csv');
disp(T1);
figure(1)
loglog(errorL2);
grid on; 
grid minor; 
end 

% function to compute the L2 for an input of u, n and f source function
function errL2 = compute_error(u,n,f)

%solve the Poisson problem first. 
[xgrid,h]=generate1d_uniform(n);
uh=fem1d(n,f);
% exact solution 
uexact=u(xgrid); 

% compute the error 
errL2 = (h^1/2)*norm(uexact(:)-uh(:),2); 
end