Need help calculating delta term in neural network back propagation

I'm trying to implement a 3 layer neural network with the following dimensions: 400 features, 40 nodes, 40 nodes, 10 targets. So, my three parameter vectors are defined as follow:

theta1 = np.random.uniform(low=0.00001, high=0.0001, size=(40,401))
theta2 = np.random.uniform(low=0.00001, high=0.0001, size=(40,41))
theta3 = np.random.uniform(low=0.00001, high=0.0001, size=(10,41))

I had no problems running forward propagation. I ran into a problem when implementing backward propagation. I get a ValueError when using np.matmul to calculate part of delta2. The code is below:

def backward_prop3(y_vectors, a1, a2, a3, a4, theta1, theta2, theta3, lamb):
   #backward propagation
   #y_vectors is vectors of results (number of rows = number of targets, number of columns = number of training examples)
   #a1 is the features, a2 and a3 are the nodes, a4 is the output
   #lambda is the
   #outputs gradient arrays for theta1 and theta2 and theta2

   m = y_vectors.shape[1]

   delta4 = a4 - y_vectors

   delta3_matmul_term = np.matmul(np.transpose(theta3),delta4)
   delta3_dot_term = np.multiply(a3, np.ones(a3.shape)-a3)
   delta3 = np.multiply(delta3_matmul_term,delta3_dot_term)
   triangle3 = np.matmul(delta4, np.transpose(a3))

   delta2_matmul_term = np.matmul(np.transpose(theta2),delta3) #problem occurred on this line
   delta2_dot_term = np.multiply(a2, np.ones(a2.shape)-a2)
   delta2 = np.multiply(delta2_matmul_term,delta2_dot_term)
   triangle2 = np.matmul(delta3, np.transpose(a2))

   triangle1 = np.matmul(delta2[1:,:],np.transpose(a1))

   
   grad3 = (1/m)*triangle3 

   grad2 = (1/m)*triangle2 

   grad1 = (1/m)*triangle1 

   return grad1, grad2, grad3

The weird thing is I had no problems calculating delta3, but ran into a problem calculating delta2. Since the code for delta 3 and delta 2 are essentially the same thing, I don't see why it worked for delta 3 and not for delta 2. Could someone suggest a possible reason please?

Edit: I have tried other numbers for the nodes, and it still doesn't work. I've also reviewed my notes, and it seem's I'm coding exactly the given equations, so I don't know what it the problem here.