I agree that the time gain from my first answer is too small to be of much help in your case. I have tested an alternative approach that may achieve time gains by a factor of 100 or more, provided that your problem satisfies the following criteria:
1) The vector p stays the same for all runs.
2) p is monotonically increasing
3) You can tolerate that the sum includes one p-value too much or too little in a small fraction of cases.
The idea is to round the limit values from dataset2 to a given number of decimal places, and convert to integers by multiplying by the appropriate power of 10. Then create an array p_index giving the index in p corresponding to all possible limit values. Although creating this array is somewhat time-consuming, I need create is only once. The same goes for array p_cum, containing the cumulative sum of p.
An example of how to generate p_index and p_um is attached as Nielsen_preparations.m.
The round off errors inevitably result in occasional errors in a small fraction of cases. This can be minimized by using more decimal places, at the cost of increasing the time for the preparations.
I can now find the relevant upper and lower indices used in the sums by direct lookup in p_index. I reduce execution time further by subtracting two values in p_cum, rather than summing a range of values.
The attached file Nielsen_test.m shows a comparison of the two methods.
