Neural network not passing gradient check

I'm following a course by Andrew ng in coursera, and while doing the course I'm also trying to implement a more sophisticated neural network along the way, so I can use what I've been learning, my neural network is a 784->10->10 2 layer network with dropout(I disabled dropout for gradient checking), but it is not passing my gradient checking approach, I've looked at it multiple times but can't find the error in the math

Result when running with gradient checking:

There is a mistake in the backward propagation! difference =
0.33333333308031154
There is a mistake in the backward propagation! difference = 0.3333333395038698
There is a mistake in the backward propagation! difference = 0.33333333720553304
There is a mistake in the backward propagation! difference = 0.3333333605889935
There is a mistake in the backward propagation! difference = 0.33333337487594783
There is a mistake in the backward propagation! difference = 0.3333333772956705
There is a mistake in the backward propagation! difference = 0.33333334149384153


#function to apply gradient checking algorithm to assert that the gradient are correct
count = 0
# Preparing necessary analytic gradient values (only require dW and dB)
for i in gradients.keys(): #Adding every weight and bias (activations and outputs excluded)
if i.startswith("dw") or i.startswith("db"):
if count == 0:
else:
count = count + 1

# Building the numerical gradient approximations
for i in parameters.keys():
for idx in np.ndindex(parameters[i].shape):
thetaplus = parameters[i][idx] + epsilon #calculating theta plus for each parameter
modified_params = parameters.copy()
modified_params[i][idx] = thetaplus
cache = self.forward_propagate(X, modified_params, self.S, True) #testing network based on modified params
J_Plus = MathUtils.cross_entropy(cache["A" + str(len(self.S) - 1)], Y)

thetaminus = parameters[i][idx] - epsilon
modified_params = parameters.copy()
modified_params[i][idx] = thetaminus
cache = self.forward_propagate(X, modified_params, self.S, True)
J_Minus = MathUtils.cross_entropy(cache["A" + str(len(self.S) - 1)], Y)
#Adding the approximation to a list
grad_approx.append((J_Plus - J_Minus) / (2 * epsilon))

#Comparing values for debugging

# Calculating relative error
difference = numerator / denominator  # Step 3'

if difference > 2e-7:
print("\033[93m" + "There is a mistake in the backward propagation! difference = " + str(
difference) + "\033[0m")
else:
print("\033[92m" + "Backward propagation works perfectly fine! difference = " + str(
difference) + "\033[0m")


Back propagation function:

    # Back propagation using Gradient descent
def back_propagate(self, X, Y, cache, parameters, S):
M = Y.shape[1] #Number of training examples
#Gradients for activations and before applying activations
gradients["dz" + str(len(S) - 1)] = (cache["A" + str(len(S) - 1)] - Y) / M
for i in range(2, len(S)):
gradients["da" + str(len(S) - i)] = np.dot(parameters["W" + str(len(S) - i + 1)].T,
gradients["dz" + str(len(S) - i + 1)])
if self.enable_dropout:
self.dropout_bw(cache, gradients, str(len(S) - i), 0.5) #Dropping out 50% of the neurons
gradients["dz" + str(len(S) - i)] = gradients["da" + str(len(S) - i)] * MathUtils.relu_deriv(
cache["Z" + str(len(S) - i)])

gradients["dw1"] = np.dot(gradients["dz1"], X.T)  # dot devido a ser a soma remember my dude produto escalar
for i in range(2, len(S)):
gradients["dw" + str(i)] = np.dot(gradients["dz" + str(i)], cache["A" + str(i - 1)].T)



Forward propagation:

    # Function to do a full forward propagation with the current parameters
def forward_propagate(self, X, params, S, train_mode):
cache = {}
#First Layer
cache["Z1"] = np.dot(params["W1"], X) + params["B1"]
cache["A1"] = MathUtils.relu(cache["Z1"])
if train_mode and self.enable_dropout:
self.dropout_fwd(cache, "1", 0.5)

for i in range(2, len(S) - 1):
cache["Z" + str(i)] = np.dot(params["W" + str(i)],
cache["A" + str(i - 1)]) + params["B" + str(i)]
cache["A" + str(i)] = MathUtils.relu(cache["Z" + str(i)])
if train_mode and self.enable_dropout: #Dropping 50% of neurons if in training mode and dropout mode
self.dropout_fwd(cache, str(i), 0.5)

# Output layer
cache["Z" + str(len(S) - 1)] = np.dot(params["W" + str(len(S) - 1)],
cache["A" + str(len(S) - 2)]) + params["B" + str(len(S) - 1)]
cache["A" + str(len(S) - 1)] = MathUtils.softmax(cache["Z" + str(len(S) - 1)])
return cache


And finally the cost function and other math functions

    @staticmethod
def cross_entropy(A, Y):
M = A.shape[1]
logprobs = np.multiply(np.log(A), Y)
cost = - np.sum(logprobs) / M
return float(np.squeeze(cost))
@staticmethod
def sigmoid(x):
return 1 / (1 + np.exp(-x))

@staticmethod
def sigmoid_deriv(x):
return MathUtils.sigmoid(x) * (1 - MathUtils.sigmoid(x))

@staticmethod
def relu(x):
return np.maximum(0, x)

@staticmethod
def relu_deriv(x):
# return np.greater(x, 0).astype(int)
x[x <= 0] = 0
x[x > 0] = 1
return x