Hexagonal Modulator Comparison with Square Modulator
For a rectangular QAM constellation, calculation of the symbol error rate is straightforward.
Exact calculation of the bit error rate is somewhat more complicated, but approximations
are more easily calculated based on the following assumptions.


Each signal point has four nearest neighbors (interior point)
or three nearest neighbors (edge points), or two nearest neighbors (corner points),
Note that 16-QAM has 4 interior points, 8 edge points, and 4 corner points.
The probability that a transmitted signal point is demodulated into 
a nearest neighbor is orders of magnitude than the probability that it is
demodulated into a more distant point
A syymbol error results in exactly one bit error (that is, data bits
are assigned to signal points in Gray code order) so that the data bits assigned
to a signal point differ from the data bits assigned to a nearest neighbor
in just one bit position.


For the 16-ary hexagonal constellation with regular equiangular hexagons
(similar to but not quite the same as the constellation shown in the question), 
the 4 interior points have 6 nearest
neighbors instead of 4 while the edge points have 4 nearest neighbors on average (some 
have 3 and some have 5, others have 4 exactly), and two corner points have 
3 nearest neighbors.  This has two effects. First, for a given spacing
between signal points, a symbol error is more likely since there are more
neighbors. Second, a symbol error can result in more than one bit error.
Thus, it is necessary to adjust the spacing between nearest neighbors
to compensate for this if one wants to have the same BER, and this will
typically increase the signal power needed. 
In short, it is not immediately obvious that


  It is supposed to give better BER measurement than Square as it has 
  wider decision boundary


is a believable claim.  In any case, the width of the decision
boundary is not quite as relevant to the BER calculation as the
distance of the constellation points from the decision boundaries,
and the number of nearest neighbors that each signal point has.



I don't decipher MATLAB very well, but if the picture with the
16 blobs was created by the MATLAB program included, then I 
don't think that the constellation that has been created is
a (regular) hexagonal constellation at all. It looks for all
the world like you began with a rectangular 16-QAM signal constellation
(which has points at $(\pm m, \pm n)$ for $m, n \in \{1, 3\}$)
and moved the second and fourth rows by $1$ to the right, so
that the $x$ coordinate is $-2, 0, 2, 4$ instead of $-3, -1, 1, 3$.
Note that this increases the maximum signal amplitude in
one branch (I in this case) from $3$ to $4$ so that 
the amplifiers have to have greater dynamic range and
the energy per bit also has increased. Any comparison
of the BER performance must take into account the fact
that the hexagonal constellation shown is operating at
a larger SNR than the standard rectangular 16-QAM constellation.
In a regular hexagonal constellation (which is what my remarks
above apply to), the points are at the
vertices (and centers) of regular (equiangular) hexagons and
are equidistant from their nearest neighbors, while
your distances are $2$ and $\sqrt{5}$ respectively.