What size n are you interested in expressing as a sum of 4 squares? The algorithm that I gave below works reasonably efficiently (say, <= 16 milliseconds) for n of many hundreds of bits. It may get bogged down if n has two or more large, distinct prime factors on the order of, say, a hundred bits. This slowness is just due to finding the prime factorization of n. If you need to work with such n, I'll implement the Rabin & Shallit algorithm (from the paper mentioned below), which doesn't use the prime factorization of n and has average time complexity O(log(n)^2*log(log(n))).

Yes, it's easy. What do you want to be the time parameter of the animation--perhaps the standard deviation? You'd likely get better pedagogy---if that's your purpose---from an Explore command than from a dedicated animation.

@vv Yes, that's the paper that I was looking for. Thank you for posting it.

When I first read this Post, I quickly suspected that the order of the permutation would be a fairly easy closed-form function of its length and that the number of disjoint cycles would be small regardless of  length; however, neither of those seem to be true. Given that, it would seem odd that Order(P(n)) <= n for all n. Has someone proven that? Below, I've verified it up to n=4096.

Here are some explorations:

Download ArabicCipher.mw

I suppose the problem could be done in Maple. Matlab's fminsearch seems equivalent to parts of Maple's third-party DirectSearch package. I need details of "Pareto front" to proceed.

@Masooma I am happy to help further if you is still need it. Newton's method and its variations are an area of expertise of mine. I still don't see any merit to your method x-> x - 2*f(x)/D(f)(x) whether for ordinary or repeated roots.

It wouldn't let me attach the worksheet above. Here it is:
 

The histograms can be generated in the worksheet. There's some trouble with uploading them directly.

Download UniformRandPrime.mw

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.