Found solution
quantiles = linspace(0,1,datasize);
[f,x, flo, fup] = ecdf(datavector); %experimental cdf
finv = interp1(f,x,quantile,'next');%Inverse experimental cdf
pdnorm = makedist('Normal'); transf = icdf(pdnorm,f); %create equivalent cdf values for N(0,1) distribution
plot(x,f); hold on; plot(finv,quantile, 'ro');