Do

solve({f(x) = y, x < x__1}, {x}, parametric, real);

This is currently only implemented for polynomial f. What you get back from that solve is a deeply nested piecewise expression which you'll need to deconstruct. For the present case, that can be done by

op([2,1,1,2], %);

It's very easy.

1. Remove from your code the line t:= 0.

2. Change the final plot command to

animate(
   fieldplot3d, 
   [[E0[x], E0[y], E0[z]], x = -.3 .. .3, y = -.3 .. .3, z = -10 .. 10, 
    grid = [10, 10, 10], arrows = SLIM, color = black
   ], t= 0..10
);

(The arguments in the outermost pair of square brackets are exactly the same arguments as your original fieldplot3d command.).

3. Re-execute the entire worksheet, including the restart.

This is just to answer your secondary question about the gridlines. I haven't yet figured out how to solve your main problem.

You can make the gridlines invisible by including option style= patchnogrid. The grid will still be in effect computationally, and the option grid= [m,n] will be effective; you just won't see the gridlines.

However, as you've already noticed, enlarging the grid won't solve your problem! It'll just make the computation lengthier and the worksheet larger and slower to load, scroll, and autosave.

Here is the code translated into Maple. I strove to maintain the character of the code (including all the variable names) while avoiding the serious shortcomings of Fortran. Despite what MacDude said, this code is not at all slow. Indeed, it's 5 times faster than RandomTools:-Generate and 6 times faster than rand(0. .. 1.), although it's 3 times slower than Statistics:-Sample.

ran1:= module()
export idum, iv, iy:= 0;
local
   ModuleLoad:= proc()
   option hfloat;
   local PARAMETER;
      Digits:= trunc(evalhf(Digits));
      PARAMETER:= 
         [
            [AM= 1./IM, NDIV= 1+iquo(IM-1, NTAB), RNMX= evalhf(1-DBL_EPSILON)],
            [IA= 16807, IM= 2147483647, NTAB= 32]
         ]
      ;
      (ModuleApply, Init):= subs(PARAMETER[], eval~([ModuleApply, Init]))[]
   end proc,

   Init:= proc(seed)
   local j;
      iv:= Vector(NTAB, datatype= integer[4]);
      idum:= max(-seed, 1);
      for j from NTAB+8 to 1 by -1 do
         idum:= irem(IA*idum, IM);
         if j <= NTAB then iv[j]:= idum end if
      end do;
      iy:= iv[1];
   end proc,

   ModuleApply:= proc(seed::integer:= 0)
   option hfloat;
   local j;
      Digits:= trunc(evalhf(Digits));
      if seed < 0 or iy = 0 then Init(seed) end if;
      idum:= irem(IA*idum, IM);
      j:= 1 + iquo(iy, NDIV);
      iy:= iv[j];
      iv[j]:= idum;
      min(AM*iy, RNMX)
   end proc
;
   ModuleLoad()
end module:

To use it, first call it with a negative integer argument to initialize it, for example, ran1(-42). After that, no argument is needed, and any nonnegative integer argument is ignored. Initializing is optional; simply starting with ran1() is equivalent to ran1(-1). A histogram of a sample of 10000 shows a reasonably uniform distribution. This can be done like this:

N:= 10000:
ran1(-42): A:= Vector(['ran1()' $ N], datatype= float[8]):
Statistics:-Histogram(A, frequencyscale= absolute);

For comparison, the other methods that I used are

B:= Vector(
   RandomTools:-Generate(list(float(range= 0. .. 1., method= uniform), N)),
   datatype= float[8]
):
C:= Statistics:-Sample(Uniform(0,1), N):
R:= rand(0. .. 1.); E:= Vector(['R()' $ N], datatype= float[8]);

As you no doubt know, the left side of your identity is Zeta(s). Since Maple also knows this, I'll use Zeta(s) for it for the rest of this article.

Constraints can be placed on variables by using an assuming clause, as in

is(Zeta(s) = product(ithprime(j)^s/(ithprime(j)^s - 1), j= 1..infinity)) assuming s > 1;

However, Maple doesn't have enough symbolic knowledge of ithprime to be able to solve this.

You can see all the symbolic knowledge (such as identities) that Maple knows about Zeta by doing

FunctionAdvisor(Zeta(s));

As far as I know, Maple has no symbolic knowledge about ithprime.

We can explore your identity through numeric computation like this (I've taken the logarithm of both sides of the identity):

nP:= 999: #number of primes to test
Ps:= hfarray([seq(ithprime(j), j= 1..nP)]):
F:= s-> evalhf(add(s*ln(Ps[k])-ln(Ps[k]^s-1), k= 1..nP)):
plot(ln@Zeta - F, 2..9);

The plot shows that the difference of the two sides of the identity is close to 0 when the infinite product is truncated.

A wrote a module to do this. It may be possible to improve the speed on this by using compiled code or by using evalhf differently. (That's not to say that this code is particularly slow!)

Pi_BaileyBorweinPlouffe:= module()
#Prints hexadecimal digits of Pi starting at an arbitrary position
local
   ModuleApply:= proc(d::posint)
   #d is the starting digit position. 9 digits will be printed because
   #that's the limit of hardware-float arithmetic according to the paper.
   option hfloat;
   local x, s, k, dm1:= d-1;
      Digits:= trunc(evalhf(Digits));
      x:= Frac(4*SSj(dm1,1) - 2*SSj(dm1,4) - SSj(dm1,5) - SSj(dm1,6));
      Digits:= 36;
      s:= sprintf("%38.36f", convert(x, binary));
      cat(seq(Hex(s[4*k-1..4*k+2]), k= 1..9))
   end proc,
      
   SSj:= proc(d::nonnegint, j::identical(1,4,5,6))
   option hfloat;
   local k;
      Frac(
         # *1
         Frac(`if`(d=0, 1./j, add((16 &^ (d-k) mod (8*k+j))/evalf(8*k+j), k= 0..d))) 
         + 
         evalhf(add(16^(d-k)/(8*k+j), k= d+1..d+14))
      )
   end proc,

   Hex:= proc(x::string)
   #Converts a string of exactly 4 binary digits to a hexadecimal digit
   option remember;
      convert(convert(x, decimal, 2), hexadecimal)
   end proc,

   Frac:= proc(x::numeric)  # *2
   option hfloat;
   local f:= x - trunc(x);
      `if`(f >= 0, f, 1+f)
   end proc
;
end module:

(* Footnotes:

*1: d=0 needs to be treated as a special case due to a bug in `mod/&^` such that
    16 &^ 0 mod m 
 returns an error. 

*2: Maple's endogenous `frac` won't work here because it may return negative. 
*)

Example usage:

Digits:= 15:
Pi_BaileyBorweinPlouffe(1);
                           243F6A888
evalf(convert(%, decimal, 16)/16^9);
                       0.141592653584667
frac(evalf(Pi));
                        0.14159265358979
Pi_BaileyBorweinPlouffe(1536);
                           362FB1341

The system and ssystem commands do the same thing as !, but they give you more flexibility for situations like this.

system(
   "curl -m 3 -o /Users/John/Documents/P\\ Admin/a.dat\"
   " --insecure https://www.google.com/?gws_rd=ssl"
);

Note carefully how I split the string over lines.

You had a semicolon at the end. I don't know whether this is required by MacOS or if you included it because you thought that it was required by Maple. If MacOS requires it, then put it at the end of the string above.

If you need the result of this command returned to Maple (as a string), then use ssystem instead.

Here's an example of generating random integers in a do loop and accumulating their lowest common multiple:

LCM:= 1:
to 9 do
   LCM:= ilcm(LCM, rand())
end do:
LCM;

The lowest common multiple is an associative operation. The technique above will work for any associative operation, and is commonly used for sums. The initial value (in this case, the 1) must be the identity element for the operation.

Having a variable, a name, that is assumed to be integer is not the same as having a value that is actually an integer (i.e., has type integer). The command numtheory:-bigomega wants an actual integer.

Now let's learn about premature evaluation, which is the essential difference between

sum(numtheory[bigomega](k), k= 2..9)

and

add(numtheory[bigomega](k), k= 2..9)

The command sum is like the vast, vast majority of other Maple commands: It evaluates its arguments before they are passed. That means that the above sum command attempts to evaluate numtheory[bigomega](k) before k has been assigned any numeric value. That causes an error.

On the other hand, the commands seq, add, and mul have special evaluation rules: Their index variable is assigned values from the prescribed range before the first argument is evaluated.

To some extent, premature evaluation can be controlled with unevaluation quotes; however, those are often difficult or impossible to use. In the above case, sum('numtheory[bigomega](k)', k= 2..9) works correctly, but it's better to just use add.

When you get an error message that says that a certain argument is expected to be a certain type, that always means that that argument must actually be that type; that error can never be corrected by making assumptions.

Your definition of N is correct. It's obvious that the midrange is approximately 10 and that the supposed answer is way off.

To get the min and max of a 2D plot, first assign the plot to a variable:

P:= plot(N, 0..20);

then

(min,max)(op([1,1], P)[.., 2]);

      9.49146814577784781, 11.3050481854892997

You can simply do

1 +~ ListTools:-PartialSums([1, 3, 4, 50, 10]);

You can't use the same variable, i, as both the index for the for loop and the index for the sum command. Simply use a different variable.

There are several other improvements that you could make, which I'll detail in a Reply.

Do

Matrix((2,2), (i,j)-> eval(p[j], solution[i]));

a) If you have actual data points, you should use Statistics:-PowerFit. It'll do the log transformation and find the closest exponent for you.

b) If R = a*P^(2/3), then R/P = a*P^(-1/3). From there, the two limits are trivial.

c) You have a lot of questions. Are you trying to do your homework for the entire semester in one night? You should indicate/show that you understand the answers that I already gave before proceeding.

Like this

restart:
K:= 100:  N__0:= 20:
NP:= R-> rsolve({N(t+1) = R*N(t)/(1+(R-1/K)*N(t)), N(0)= N__0}, N(t), makeproc):
plot(
   [seq(['[k, NP(R)(k)]' $ k= 1..10], R= [2, 5, 10])],
   legend= (R=~ [2, 5, 10]),
   labels= [t,N]
);