Does anybody know where to start modeling a square multimode fiber with top hat output?

I have attached a paper from Optoskand where they describe a fiber with a square core. Where the beam exits the fiber the intensity profile has a square top hat distribution.

The beam has diverges out of the fiber, where the divergence angle is given by the BPP (beam parameter product) and core size.

The beam propagating out of the fiber loses its top hat distribution. Collimating this beam and focusing it back in to a spot results again in a top hat distribution in the focal plane.

What would be a good approach to model such a beam?