Thank you so much for the detailed responses. 
So, my approach to calculate relu’s gradient was flawed and an alternative but inefficient way of calculating grad_h would be:
    grad_h = grad_h_relu.copy()
    grad_h = grad_h_relu.dot(Jacobian(h_relu))

right?
Also, in this example, we are processing a batch of 64 examples together. h is of size [64, 100] and so is h_relu. So how do we calculate the Jacobian of matrices?