Unexpected result in finite difference solution beam deflection problem

Question

I try to numerically solve (by finite difference method) the beam deflection equation $$EJ\cdot\frac{\partial^4 v}{\partial x^4} = -g\cdot F\cdot\rho$$ with following boundary conditions $$v(0)=0;\frac{\partial v(0)}{\partial x} = 0$$ $$\frac{\partial^2 v(L)}{\partial x^2} = 0;\frac{\partial^3 v(L)}{\partial x^3} = 0$$

I'm using following formulas to approximate equation and its boundary conditions: $$\frac{\partial v}{\partial x} = \frac{v_{i+1} - v_{i-1}}{2h}$$ $$\frac{\partial^2 v}{\partial x^2} = \frac{v_{i+1} -2v_i + v_{i-1}}{h^2}$$ $$\frac{\partial^3 v}{\partial x^3} = \frac{-v_{i-2} + 2v_{i-1} - 2v_{i+1} + v_{i+2}}{2h^3}$$ $$\frac{\partial^4 v}{\partial x^4} = \frac{v_{i-2} - 4v_{i-1} + 6v_i - 4v_{i+1} + v_{i+2}}{h^4}$$

For check result I use the simple analitical formula $$\delta_{max} = \frac{\omega \cdot L^4}{8EJ}$$

Here is my Python code for it problem

import numpy as np
import matplotlib.pyplot as plt


E = 201e9
rho = 7850
G = 9.81
d = 0.1
F = 0.25 * np.pi * d ** 2
J = np.pi * d ** 4 / 64
L = 3
N = 100
#N = 1000
#N = 5000
#N = 7000

A = np.zeros((N+1, N+1))
f = np.zeros(N+1)
mesh, dx = np.linspace(0, L, N+1, retstep=True)


omega = -G * rho * F

delta_max = omega * L ** 4 / (8 * E * J)
theta_end = omega * L ** 3 / (6 * E * J)

print(delta_max * 1e3)



for i in range(2, A.shape[0]-2):
    A[i, i] = 6 * E * J
    A[i, i-1] = -4 * E * J
    A[i, i+1] = -4 * E * J
    A[i, i-2] = E * J
    A[i, i+2] = E * J
    f[i] = -G * rho * F * dx ** 4


A[0, 0] = 1
f[0] = 0

A[1, 0] = -4 * E * J
A[1, 1] = 7 * E * J
A[1, 2] = -4 * E * J
A[1, 3] = E * J
f[1] = -G * rho * F * dx ** 4

A[-2, -1] = -2 * E * J
A[-2, -2] = 5 * E * J
A[-2, -3] = -4 * E * J
A[-2, -4] = E * J
f[-2] = -G * rho * F * dx ** 4

A[-1, -1] = 2 * E * J
A[-1, -2] = -4 * E * J
A[-1, -3] = 2 * E * J
f[-1] = -G * rho * F * dx ** 4


v = np.linalg.solve(A, f)


fig, ax = plt.subplots()
ax.plot(mesh, v*1e3)
plt.show()

delta_max_calc = v[-1]

theta_end_calc = np.gradient(v, dx, edge_order=2)[-1]

print(100 * (delta_max - delta_max_calc) / delta_max)

print(100 * (theta_end - theta_end_calc) / theta_end)

The problem is: the error increase when i increase number of cells in my simulation.

I get following results: $err_{100 \space cells} = -0.009\%$, $err_{1000 \space cells} = 0.0038\%$, $err_{5000 \space cells} = 2.34\%$, $err_{7000 \space cells} = 8.86\%$

I'm not very good in numerical analysis, so the question is: is this result is correct or i incorrectly approximate boundary conditions?