The error in your procedure is that ExitCond needs to be reinitialized to false for each value of k.

A minor issue is that a procedure should declare the modules (or packages) that it uses. Make the first line of the procedure

uses LinearAlgebra, Statistics;

Here is Maple code for the above algorithm. The return value is a list of the m-1 centers (each a Vector), followed by the common radius. You did not say that the balls needed to lie entirely within the bounds, merely that their centers did, so I have not required that.

Balls:= proc(
   Box::range(Vector(realcons)), #range for each dimension
   X1::Vector(realcons), #initial ball's center
   R1::And(realcons, positive), #initial ball's radius
   m::posint #number of balls (incl. initial)
)
uses LA= LinearAlgebra, C= combinat;
local
   d:= numelems(X1), #number of dimensions
   L:= lhs(Box), #lower bounds
   U:= rhs(Box), #upper bounds
   M:= U-L, #magnitudes of dimensions 
   RV:= proc()
   local 
      RV:= (()-> LA:-RandomVector(d, 'generator'= 0. .. 1.) *~ M + L),
      rv:= RV()
   ;
      while LA:-Norm(rv - X1, 2) <= R1 do rv:= RV() end do;
      rv
   end proc,
   X:= ['RV'() $ m-1]
;
   X, 
   min(
      LA:-Norm~((`-`@op)~(C:-choose(X,2)), 2)[], #pairwise distances
      (LA:-Norm~(map(`-`, X, X1), 2) -~ R1)[] #distances from initial ball
   ) / 3
end proc
:

Example usage:

Balls(<-1,-1,-1>..<1,1,1>, <.1,.2,.3>, 0.6, 15);

Would someone please define alpha and alpha channel for me?

@Carl Love The Question has been changed a few times since I posted my Answer above. My Answer pertained to something like

plot(-x^2+1, x= -1..1);

which was at the end of the Question as an example of something much simpler that gave the same error message. It was this that I was responding to. This line has no error.

@Jayaprakash J No, you do not need Mathematica to do this.

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?