Yes, you are right and I have assumed you are passing an already reset bn layer into your custom implementation.
Your MyBatchNorm2d implementation seems to work correctly:
class MyBatchNorm2d(nn.Module):
    def __init__(self, init_bn):
        super(MyBatchNorm2d, self).__init__()
        self.bn = copy.deepcopy(init_bn)
        
        # Aux BN
        self.bn_aux = copy.deepcopy(init_bn)
        self.bn_aux.reset_parameters()
        
        self.register_buffer('num_batches_tracked', torch.tensor(0, dtype=torch.long))
        
    def forward(self, input):            
        if self.training:
            self.num_batches_tracked = self.num_batches_tracked + 1            
        
        # Odd runs are AUX input
        if self.training and self.num_batches_tracked % 2 != 0:
            return self.bn_aux(input)
      
        return self.bn(input)


bn = nn.BatchNorm2d(3)

with torch.no_grad():
    bn.weight.add_(10)
    bn.bias.add_(10)
    bn.running_mean.add_(10)
    bn.running_var.add_(10)

my_bn = MyBatchNorm2d(copy.deepcopy(bn))

x = torch.randn(2, 3, 4, 4)
out1 = bn(x)
out2 = my_bn(x)
print(out1)
print(out2)

How are you comparing the results using the test error?