Why is my laser beam intensity profile not Gaussian?
What could cause such a shape?
  Could it be due to the laser source which could not be perfectly single-mode?
  Or could it be due to the sensor?
  


There is nothing particularly wrong with the "shape", but there are a few things you can do on the sensor and data processing side, to improve the extraction of an accurate profile.

Your "biggest" problem is going to be the saturation of your web cam's sensor because of the intensity of the beam. This is evident as a flattening of the Gaussian curve's peak. There is a bit of that in this example, but if the $Y$ axis is not scaled, then values in the vicinity of 50 don't look like being saturated....However...

Please make sure that "Automatic Exposure Control" on the web cam is off and then manually set the exposure for minimum distortion (it will probably have to be set to its lowest setting). This will also affect the slopes of the curve as exposure control will strive to maintain contrast across the image.

In addition, you might still need some "analogue" attenuation too. To achieve this, use one or more neutral density filters in front of the camera.

Finally, as it is evident from this profile, it is composed of integer quantities and you are trying to fit a continuous distribution on it. There are two things that can be done here towards improvement: The first, is to shoot a video and then create an average image of the profile over a few hundred frames. This will give you a much smoother real ($\mathbb{R}$) number profile which will vastly improve the fit. The second, is to still shoot a video but produce the beam profile as a boxplot. This will give you the opportunity to observe exactly what the profile "does" over a few hundred (or more) frames. Very briefly, a boxplot would present in one view, information about the distribution of the brightness values at each pixel of the profile over time. You can then fit a Gaussian over the means (or medians) of the profile.

Hope this helps.