How to plot the pressure distribution over an airfoil?
I have found only the description of this pattern with experimental graphs.


That's because that is the best way to have detailed and precise data: you either simulate (e.g., with Nastran-Patran or XFoil, as cited by Peter Kämpf in his answer) or build a [scale] model and use a wind tunnel. If you use a scale model, depending on the focus of the experiment, you will have to be careful (while scaling back up) not to damage your data: for this you have adymensional quantities to look after, such as the Reynolds number.


  Could anyone please tell me the mathematical formula required to plot the distribution?


There is not one formula to do that, but rather a method that given the airfoil shape (and the angle of attack) will estimate the pressure distribution. This method was extensively used in the days before computers and FEM tools existed (for what I know SR71 is a notable example where this method has been used in the early phase of the wings design), but, for what I know, nowadays is seldom used.

The method is based on Conformal Mapping, in particular the Joukowsky transform or the Kármán–Trefftz for more detailed and general cases (see Sec. 4.2 of this pdf for a more in-depth explanation). This is a transformation in the complex plane that, from the airfoil profile that you want to analyse, will map it onto a circle.

You want to do this because the ideal flow around a circle is well known:


At this point you have to know that around an airfoil the air will have a certain amount of circulation and that the total flow around it is the sum of the ideal unviscous flow plus the flow induced by the circulation:



 Image from MIT 

For this reason, you want to compute the ideal flow around a cilinder that has a rotation speed equivalent to the circulation around your airfoil[*]:



and it turns out that we have some formulas for that (image above taken from this page).



[*]: note that you will have different circulations for different angles of attack.



Now, why did we do all this? Because once we have the flow around the rotating cilinder and 
we have the mapping from the airfoil to the cylinder, we "only" have to invert the mapping:



and we immediately have our flow around the airfoil.

And once we have the airflow, finding the pressure distribution is quite straightforward, thanks to Bernoulli (you have the velocity distribution and you can assume that the transit of the air around the wing is fast enough to be adiabatic)