A small correction: `commonind` returns the first common Index found, but doesn't error if there are more than one. From the code you posted, you are setting `M[b+j] = L`, which looks right. However, it also seems like you should be setting `M[b+j+1] = R` within your loop, could that be part of the problem you are seeing?
In the original `replacebond(...)` function, the point of the line `tags = tags(linkind(M, b))` is so that the new link index that is being introduced by `factorize` has the same tags as the original link index. Likely you would want this feature in your own function, but it isn't absolutely necessary. However, if the call to `tags(linkind(M, b+j))` is failing as you say, then the function `linkind(M, b+j)` is not finding a link index between your MPS tensors `b+j` and `b+j+1`. This would indicate that there is something off with your MPS, since it should have one common Index between tensors `M[b+j]` and `M[b+j+1]`. As Miles recommended, you should print the indices of those tensors to try to debug what is going wrong.
The `setprime` lines in the `factorize` function is there so that you can set the prime level of the new index being introduced if you want. `replacebond` doesn't set this, but other applications might want to. You should think of `factorize` as a more general functionality which is used in a wide variety of applications, and happens to be used in `replacebond` as a specific application. The final return of the Index `l` from factorize is just for convenience, if you want a different way to "grab" the new Index that is introduced by factorize (it returns the Index `l == commonind(L, R)`). For the use in `replacebond` it is just being discarded.
Also, this week I will be working on adding functionality to ITensors.jl for a generalization of `replacebond` to a general number of sites. Feel free to continue working on your own version, but hopefully there will be built-in functionality for that soon!