You are right it should probably be something like,
PTDF[bus_name, branch_name] = b_val*(B_inv.matrix[bus_fr, n] - B_inv.matrix[bus_to, n])
But I did not check this in detail. My goal was to show you just roughly how it was done. My recommendation would be to take the precise definition out of a textbook to make sure.
Thanks for the tip about AxisArrays, I had missed the atvalue function!