As to why your results aren’t matching the ones in the article, does the calculation in the article use “hard core” bosons restricted to have a maximum occupancy of 1, or regular bosons which can (in principle) have any number of bosons on a site?
I ask because your code looks like it is just using our “Boson” site type without any “dim” parameter, in which case the maximum occupancy defaults to just 2. If you are trying to approach the limit of bosons with higher occupancies, it’s important to raise the maximum dimension of the boson sites until you can confirm that your results are converged with this dimension value.
You can make more general sites this way:
boson_dim = 4
sites = [isodd(j) ? siteind("Boson";dim=boson_dim) : siteind("S=1/2") for j=1:N]
If you do @show sites you will see that now all of your Boson sites will have a dimension of 4, and generally whatever maximum dimension you set.