@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.

@acer Shouldn't your `mod/ceil` procedure have ceil(a) where it now has ceil(x)?

@Rouben Rostamian  There are a few more problems:

1. The last two initial/boundary conditions aren't orthogonal to the boundary, i.e., they should each specify the other independent variable.
2. With these conditions, you need to specify the range of the space variable. One of the two endpoints of the range must be 0.
3. You may need to specify which of the independent variables is the time variable; not sure. It certainly wouldn't hurt to do so.

The diagram in the link says nothing about A and B other than that they are on that circle of intersection. More-specific information is needed.