We present an algorithm of finding numerical solutions of pulsar equation. The problem of finding the solutions was reduced to finding expansion coefficients of the source term of the equation in a base of orthogo- nal functions defined on the unit interval by minimizing a multi-variable mismatch function defined on the light cylinder. We applied the algorithm to Scharlemann & Wagoner boundary conditions by which a smooth solu- tion is reconstructed that by construction passes success- fully the Gruzinov's test of the source function exponent.

The Monte Carlo simulations of pulsar periods, pulse-widths and magnetic inclination angles are performed. Using the available observational data sets we study a possible trial parent distribution functions by means of the Kolmogorov-Smirnov significance tests. We also use an additional condition that the numbers of generated interpulses, whether from both magnetic poles or from single pole, are at the observed levels. We conclude that the parent distribution function of magnetic inclination angles is neither flat nor cosine but it is a more complicated function with a local maximum near alpha=25deg and another weaker one near alpha=90deg. The plausible distribution function of pulsar periods is represented by the gamma function. The beaming fraction describing the fraction of observable radio pulsars is about 0.12.