It would be better to totally omit the output from your Questions than to include it unformatted. The inclusion of the output makes the input more difficult to see.

I saw a few minutes ago that you had another question, about pseudo-differential operators. But that's gone now. Please repost that, preferably in a new Question thread.

@Mayo_Kun What are the units, or at least the dimensions, of f, of theta, and of eta? This is a fluid dynamics problem, right? There must be something in the physics of the problem that allows the 2nd derivative of f to be raised to the 0.4 power.

@Masooma Your modified Newton's method (x - 2*f(x)/D(f)(x)) seems to converge rarely for the cases that I've tried, and even if it converges, it does so slowly. Do you have some reference for this method? Are you sure that you read it correctly?

@vv Yes, admittedly, my first procedure EasyRandPrime is biased towards primes that are at greater distances from their predecessors. That's why I also provided UniformRandPrime.

@vv It seems like you haven't run my code!  7 and 11 are consecutive primes:

seq(UniformRandPrime(7..11), k= 1..30); 
7, 7, 7, 11, 7, 7, 7, 11, 11, 7, 7, 7, 11, 11, 11, 7, 7, 7, 11, 11, 11, 11, 11, 11, 11, 7, 11, 7, 7, 7

That looks uniform to me.

@Mohamed19 So, if your formula is correct then both Mathematica and Maple are wrong and give the same wrong formula?

@Masooma Okay, I now understand what you mean by multiple roots; I call them repeated roots. Can you provide a reference that shows that the method that you propose has order of convergence 2 for repeated roots?

An effective way to find repeated roots is to divide the function by its derivative, which converts them to simple roots. Then apply any method that converges well for simple roots.

@tomleslie I suspect that you're already well aware of this, but it's worth pointing out for the benefit of the OP and other readers: The method that you propose is only feasible when the upper limit of the primes selected from is reasonably small, say less than a few million. In those cases, your method is indeed reasonably efficient and perfectly uniform. However, some applications, crytography in particular, require much much larger random primes. In those cases, your method is totally infeasible.

@Mohamed19 Which do you mean by "this derivative"? Do you mean that you want to define a function whose derivative is the incorrect formula that you gave? That would be possible, but is it worth it?

@Carl Love The Wolfram (essentially, Mathematica) reference page for BesselI gives three (presumably equivalent) formulations of diff(BesselI(nu, z), z). One of those is the same as what Maple gives.

@Mohamed19 Let's suppose that your derivative were correct. That would imply that the function is not differentiable for alpha = -1. But the plot

plot(BesselI(-1, x), x= -2..2);

shows that it is differentiable.

@Mohamed19 Consider the case alpha = -1. If your expression for the derivative were correct, it would imply that BesselI(-1, x) was a constant or piecewise-constant function. Simple plotting shows that that's not true.

Maple gives the derivative (wrt x) as

diff(BesselI(alpha+1, x), x);
BesselI(alpha+1,x)+alpha/x*BesselI(alpha,x)
 

@Mohamed19 

The derivative of order n can be computed like this:

diff(BesselI(alpha, x), [x$n]);

This works for any nonnegative integer n or for unspecfied n.

An infinite series representation of the form x^alpha*P(x) with P(x) a power series can be obtained by 

convert(BesselI(alpha, x), Sum);

@Adam Ledger I don't see your point. The goal is to get a square wave, not a staircase.