Sol:= dsolve({diff(y(x),x)=y(x), y(0)=1}, numeric, output= listprocedure):
theta:= eval(y(x), Sol):
evalf(Int(ln@theta, 0..2));
                   1.999999704723899

The exact answer is 2. So the accuracy is limited to the accuracy of dsolve not the accuracy of numerical integration.

Is this what you had in mind for the second plot?

restart:
b:= 100:
plot([seq(BesselJ(0, BesselJZeros(0,n)*sqrt(z/b)), n= 1..4)], z= 0..100);

And for the animation, I had to assume some values for your constants.

restart:
A:= K[n]*cos(BesselJZeros(0,n)^2/4/b*sqrt(g)*t - sigma[n])*BesselJ(0, BesselJZeros(0,n)*sqrt(z/b)):
n:= 3: K[n]:= 1: b:= 100: g:= 1: sigma[3]:= Pi/3:
plots:-animate(plot, [A, z= 0..100], t= 0..20);

The closest thing to that that is possible is

(A,C,B,F):= convert(Q, list)[];

Note that the order is not (A,B,C,F).

Your list6 is too big to generate the permutations all at once, which is what combinat:-permute does. You need to use the iterator commands, which generate them one at a time. These are combinat:-firstperm, combinat:-nextperm, etc.

You must mean for the cumulative probability to be greater than 0.95, not the one-point probability. For the cumulative probability, you need the sum of the one-point probabilities from k = 0 to n.

restart:
lambda:= 2:
sum(exp(-lambda)*lambda^k/k!, k= 0..n):
fsolve(% = .95, n);
                        4.04766491450491

Since n must be integer, we take n = 5.

Download dchange.mw

That's one solution. There are several others. Obviously < 0, 0, 0 > is a solution.

Download MultiNewton.mw

You have two separate problems. The problem that is preventing the plotting is that Pi in Maple is spelled with a capital P when you mean the well-known mathematical constant; lowercase pi is just a variable that prints the same way.

The problem with s(2) is that you haven't made s a procedure; rather, s is an expression. To make it a procedure, do

s:= theta-> piecewise(...);

and change the plot command to plot(s(theta), theta= 0.01..2*Pi) or plot(s, 0.01..2*Pi).

If you want to keep s as an expression, then to evaluate it at 2 do

eval(s, theta= 2);

Just use the collect command.

P:= (1+1/z)^2*(1+z)^2*(r[1]*z+r[0]+r[1]/z):
collect(%, z);

Maple does not ordinarily show negative exponents, using denominators instead; although the negative exponents do appear in the internal representation.

In some cases, it is necessary to use expand before collect.

Use dsolve followed by odeplot.

The next problem is your initial conditions. Clearly there is no solution for u1(0)=0, since u1(t) appears as a denominator in the system.

Sol:= dsolve(sys union {u1(0)=1, u2(0)=1, u3(0)=1}, range= 0..15, numeric):
plots:-odeplot(Sol, [u1(t),u2(t),u3(t)], t= 0..15);

A key observation is that an upper bound for the exponents is the base-2 log of the largest element.

f:= proc(ns::seq(And(posint, Not(identical(1)))), $)
local Ns:= [ns], M:= ilog2(max(Ns)), k, R:= Vector(1, [nops(Ns)]);
     if Ns=[] then  return []  end if;
     for k from 2 to M do
          R(k):= nops(select(n-> iroot(n,k)^k=n, Ns))
     end do;
     for k from M by -1 while R(k)=0 do  R(k):= ()  end do;
     convert(R, list);
end proc:

f(27);
                           [1, 0, 1]
f(2,3,4,9,81,1024);
                 [6, 4, 0, 1, 1, 0, 0, 0, 0, 1]

Did you recently study a section on Bezier curves?

Start with a piece of graph paper. Write your name in large letters as simply as possible. Designate an arbitrary point on the paper as the origin (0,0). Mark all the significant points: endpoints, intersections, sharp turns. Get their coordinates. For the straight sections, just plot as the line segment between two points. For the curved pieces, pick two more "control points" in addition to the endpoints (the control points are not necessarily on the curves). Parmetrically plot the unique cubic Bezier curve for each group of four points (for each curve). Then put the plots together with plots:-display.

Your function definition doesn't work because the syntax of a function definition is

f:= (x,y)-> ...

not

f(x,y):= ....

But, instead of a function, it would be much easier to make f an Array, like this:

LL:= [[-11, -6, -1, 5, 5, 5, 5, 5], [-6, 1, 3, 6, 9, 9, 9, 9], [0, 5, 10, 13, 8, 6, 6, 6],
[9, 11, 12, 10, 10, 9, 9, 9], [9, 10, 11, 12, 12, 5, 5, 5], [11, 11, 11, 11, 10, 7, 7, 7],
[12, 12, 11, 13, 12, 7, 7, 7], [12, 12, 11, 11, 13, 8, 8, 8]]:
f:= 150+Array(0..7, 0..7, ListTools:-Transpose(LL)):
add(f[0,y], y= 0..7);
                              1236

Note that to save typing I made each original entry the difference of the number from 150. Then I added 150 to everything at the end.

Here's how (one way) to enter this system of PDEs into Maple. Like Preben, I need some clarification on the initial conditions.

restart:
alias(X= x(t,tau), Y= y(t,tau), Z= z(t,tau), Q= q(t,tau)):
xd,yd,zd,qd:= map(diff, [X,Y,Z,Q], t)[]:
xt,yt,zt:= map(diff, [X,Y,Z], tau)[]:
M:= (xt*(-Y-Z)+yt*(X+a*Y)+zt*(b-Z*(X-mu)))/(xt^2+yt^2+zt^2):
XM:= -Y-Z-xt*M:  YM:= X+a*Y-yt*M:  ZM:= b+Z*(X-mu)-zt*M:
pde1:= xd = XM+epsilon[x]*Q:
pde2:= yd = YM+epsilon[y]*Q:
pde3:= zd = ZM+epsilon[z]*Q:
pde4:= qd = a__x*XM+a__y*YM+a__z*ZM+beta*Q:

#Boundary conditions
BCs:= x(t,0)=x(t,2*Pi), y(t,0)=y(t,2*Pi),
     z(t,0)=z(t,2*Pi), q(t,0)=q(t,2*Pi),
     D[2](x)(t,0)=D[2](x)(t,2*Pi), D[2](y)(t,0)=D[2](y)(t,2*Pi),
     D[2](z)(t,0)=D[2](z)(t,2*Pi), D[2](q)(t,0)=D[2](q)(t,2*Pi)
:
#Initial conditions
# Need clarification on these.
ICs:= NULL:

#Parameters:
mu:= 1.85:  beta:= -2.23:  a:= 1/2:  b:= 3/4:
epsilon[x]:= -1/2:  epsilon[y]:= -1:  epsilon[z]:= -1/2:
a__x:= -1:  a__y:= 2:  a__z:= -1:

#Get the solution.
#The command is pdsolve, not dsolve, since these are PDEs, not ODEs.
Sol:= pdsolve({pde||(1..4)}, {BCs, ICs}, numeric):

Error, (in pdsolve/numeric) initial/boundary conditions must be defined at one or two points for each independent variable

How about

ex:= f(h*a[3]+x, y+h(a[3]-b[32])*k[1]+h*b[32]*k[2]):
mtaylor(ex, [x=0, y=0]);