Set the differentiation matrix as A to be N+2 x N+2
A = np.zeros((N+2,N+2))
Here we set them using a loop, which must be N+1 long (the size of the off diagonals)
for i in range(N+1):
    A[i,i+1] =  1./(2*h)
    A[i+1,i] = -1./(2*h)
    
Does @ mean multiply in this line below?
df_cent = A@f(x) # do the multiplication
 
    