Context : Aims : The study of accurate methods to estimate the distribution of stellar rotational velocities is important for understanding many aspects of stellar evolution . From such observations we obtain the projected rotational speed ( v \sin i ) in order to recover the true distribution of the rotational velocity . To that end , we need to solve a difficult inverse problem that can be posed as a Fredholm integral of the first kind Methods : In this work we have used a novel approach based on Maximum likelihood ( ML ) estimation to obtain an approximation of the true rotational velocity probability density function expressed as a sum of known distribution families . In our proposal , the measurements have been treated as random variables drawn from the projected rotational velocity probability density function . We analyzed the case of Maxwellian sum approximation , where we estimated the parameters that define the sum of distributions . Results : The performance of the proposed method is analyzed using Monte Carlo simulations considering two theoretical cases for the probability density function of the true rotational stellar velocities : i ) an unimodal Maxwellian probability density distribution and ii ) a bimodal Maxwellian probability density distribution . The results show that the proposed method yielded more accurate estimates in comparison with the Tikhonov regularization method , especially for small sample length ( N = 50 ) . Our proposal was evaluated using real data from three sets of measurements , and our findings were validated using three statistical tests . Conclusions : The ML approach with Maxwellian sum approximation is a accurate method to deconvolve the rotational velocity probability density function , even when the sample length is small ( N = 50 ) .