I have first four cumulants i.e. mean, variance, skewness and kurtosis. Is there any way to find probability density function, numerically or analytically?
In MATLAB one can find normal probability distribution using first two cumulants by normpdf(x,mean,varaince) function.
On
There is no probability distribution other than the Gaussian having finite number of non-zero cumulants. This was a theorem proved by Marcinkiewicz. See http://www.scholarpedia.org/article/Cumulants#Uniqueness
So you cannot reconstruct a distribution with just the first four cumulants, not because it is not unique like in the case of Taylor expansions, but because no such distribution even exists.
If your question is: "I have a random variable whose mean is $a$, variance is $b$, skewness is $c$ and kurtosis is $d$. Can I determine the PDF?" Then the answer is no, for essentially the same reason that a function is not determined by only its first four Taylor coefficients.
If, on the other hand, you know the type of random variable but have some unknown parameters -- e.g. your variable is $\textrm{Uniform}(r, s)$ or $\textrm{Normal}(\mu, \sigma)$ -- and you know the first four cumulants, then this often is enough information to infer the missing parameters.