Back propagation on matrix of weights

I am trying to implement a Neural Network for binary classification using python and numpy only.

My network structure is as follows: input features: 2 [1X2] matrix Hidden layer1: 5 neurons [2X5] matrix Hidden layer2: 5 neurons [5X5] matrix Output layer: 1 neuron [5X1]matrix I have used the sigmoid activation function in all the layers.

Now lets say I use binary cross entropy as my loss function. How do I do the back propagation on these matrices to update weights?

class Layer():
    def __init__(self,number_of_neurons,number_of_inputs):
        self.weights=np.random.rand(number_of_neurons, number_of_inputs)
        self.bias=np.random.rand(number_of_neurons,1)     

class NeralNetwork():
        def __init__(self, layer1, layer2,layer3):
            self.layer1 = layer1
            self.layer2 = layer2
            self.layer3 = layer3

    def sigmoid(self,x):
        return 1 / (1 + np.exp(-x))

    def derivative_sigmoid(self,x):
        return x*(1-x)

    def get_cost_value(self,Y_hat, Y):
        m = Y_hat.shape[1]
        cost = -1 / m * (np.dot(Y, np.log(Y_hat).T) + np.dot(1 - Y, np.log(1 - Y_hat).T))
        return np.squeeze(cost)

    def get_cost_derivative(self,Y_hat,Y):
        return  - (np.divide(Y, Y_hat) - np.divide(1 - Y, 1 - Y_hat))

    def train(self,inputs,labels,epocs):
        for epoc in range(1,epocs+1):
            z1=np.dot(self.layer1.weights,inputs)+self.layer1.bias
            a1=self.sigmoid(z1)
            z2=np.dot(self.layer2.weights,a1)+self.layer2.bias
            a2=self.sigmoid(z2)
            #print(a2.shape)
            z3=np.dot(self.layer3.weights,a2)+self.layer3.bias
            a3=self.sigmoid(z3)
            #print(a3.shape)
            if epoc%100 is 0:
                print(a3)

            cost=self.get_cost_value(a3,labels)
            #print(cost)



            layer3_delta=self.derivative_sigmoid(a3)*self.get_cost_derivative(a3,labels)
            print(layer3_delta.shape)
            Dw_layer3=np.dot(layer3_delta,a2.T)
            Db_layer3=layer3_delta
            #print(Dw_layer3.shape)
            layer2_delta=np.dot(self.layer3.weights.T,layer3_delta)*self.derivative_sigmoid(a2)
            #print(layer2_delta.shape)
            Dw_layer2=np.dot(layer2_delta,a1.T)
            Db_layer2=layer2_delta

            layer1_delta=np.dot(self.layer2.weights.T,layer2_delta)*self.derivative_sigmoid(a1)
            Dw_layer1=np.dot(layer1_delta,inputs.T)
            Db_layer1=layer1_delta
            #print(Dw_layer1)

            self.layer1.weights-=((1/epoc)*Dw_layer1)
            self.layer2.weights-=((1/epoc)*Dw_layer2)
            self.layer3.weights-=((1/epoc)*Dw_layer3)

            self.layer1.bias-=((1/epoc)*Db_layer1)
            self.layer2.bias-=((1/epoc)*Db_layer2)
            self.layer3.bias-=((1/epoc)*Db_layer3)

So far I have tried to implement this as shown above. But there seems to be a mistake because, after training, the network doesn't seem to have learned. Please let me know if you have any inputs.