Make it a single command:

eval(y, isolate(x, D(a)));

Do you want a uniform random selection from all 6-tuples of nonnegative integers whose sum is 8? It can be done like this:

C:= combinat:-composition(8+6, 6):
C[(rand() mod nops(C))+1] -~ 1;

If you want to make such a random selection repeatedly, the above can be made more efficient. Let me know.

Your attempt to trap division-by-zero errors went awry. You set it up so that division by zero returns infinity, but the plots:-circle command doesn't know what to do with infinity as an argument. Your problem can solved by simply filtering out the values that lead to division by zero; there is no need to trap those errors. 

To fix the problem, simply replace the line

d := seq(circle([1, 1/i], 1/i, color = blue), i = -20 .. 20, .1);

with

d:= seq(`if`(i=0, NULL, circle([1, 1/i], 1/i, color= blue)), i= -20 .. 20, .1):

The principles for solving the system of piecewise ODEs are the same as the simplified example above. It is done like this:

Sys1:= EQ||(1..3), op(2, EQ4), op(2, EQ5):
Sys2:= EQ||(1..3), op(4, EQ4), op(4, EQ5):
ICs1:= seq(q||k(0)=0, k= 1..5), seq(D(q||k)(0)=0, k= 1..3):
Sol1:= dsolve({Sys1, ICs1}, numeric, maxfun=0, range= 0..10):
ICs2:= eval(
     [seq(q||k(10) = q||k(t), k= 1..5),
      seq(D(q||k)(10) = diff(q||k(t), t), k= 1..3)
     ], Sol1(10)
)[]:
Sol2:= dsolve({Sys2, ICs2}, numeric, maxfun= 0):
P1:= plots:-odeplot(Sol1, [seq([t, q||k(t)], k= 1..5)], t= 0..10):
P2:= plots:-odeplot(Sol2, [seq([t, q||k(t)], k= 1..5)], t= 10..20):
plots:-display([P1,P2]);

The solution to this ODE becomes singular before t = 0.03. However, I will show how to work with an ODE like this using the function f(x,t) from your earlier post. That one doesn't go singular until about t = 15.

Individual variables can be cleared with unassign.

This may not be the most efficient way, but I think that it is robust and that it handles all the special cases.

power:= (term,x)-> `if`(term=0, undefined, expand(diff(term,x)*x/term));

Use piecewise:

f:= t-> piecewise(t < 10, 1-t, t):

f(t);

You have not defined the function for t <= 0 or t = 10. The above will of course apply the upper branch 1-t to t <= 0 and the lower branch to t = 10.

Try this, if you haven't already:

for T from 0 to 15 by 0.1 do
     your script
end do;

If that doesn't work, then you'll need to post your code.

Yes. The default domain of computation in Maple is the complex numbers; you do not need to make any special declaration.

The complex unit---the square root of -1---is by default represented in Maple by capital I. This can be changed if needed.

The command evalc attempts to break an expression into its real and imaginary parts assuming that any variables in the expression represent real numbers, unless the variables are explicitly declared otherwise. Example:

expr:= exp(x+I*y);

evalc(expr);

evalc(expr) assuming y::complex;

Note that the assumptions that x and y are real are only active during the execution of the evalc command.

The commands Re and Im extract the real and imaginary parts, respectively, of an expression without making any assumptions. They are often used in combination with evalc. Example:

Im(expr);

evalc(Im(expr));

First, I factored out the constant 2*gamma*cos(psi)*cos(theta). Doing this does not affect Maple's ability to do the integral; it just makes the presentation neater. Then I switched the order of integration and added assumptions D > L and L > 0.

j:= (1/(D-r*sin(theta)*sin(phi))-1/(D+r*sin(theta)*sin(phi)))*r^2*sin(phi);

int(int(j, phi= 0..Pi), r= 0..L) assuming D > L, L > 0;

Doing it with

PolynomialTools:-CoefficientVector(p, x, termorder= reverse, storage= sparse)

is much more efficient than using the naive method of calling coeff iteratively.

restart:
NaiveCoeffs:= (p::polynom, x::name)->
     map[3](coeff, p, x, [seq(degree(p,x)..0,-1)]):
N:= 2^11:
P:= ['randpoly(x, degree= N)' $ N]:
time(map(NaiveCoeffs, P, x));

     2.218

time(map(PolynomialTools:-CoefficientVector, P, x, termorder= reverse, storage= sparse));

     0.109

This is a DAE system because the EQ3 has no derivatives. Since your previous question had a similar system, which was stiff, I assumed this one was stiff. So I used the stiff DAE solver method= rosenbrock_dae. This produced good, quick results.

Sol:= dsolve(
     {EQ||(1..5), seq(q||k(0)=0, k= {1,2,4,5}), D(q1)(0)=0, D(q2)(0)=0},
     numeric,
     method= rosenbrock_dae
):
plots:-odeplot(Sol, [seq([t,q||k(t)], k= 1..5)], t= 0..0.1);

You can use Maple's plot command to generate the Matrix that you export. This seems useful because you are exporting the Matrix for the purpose of plotting. Maple's plot will automatically choose sufficient intermediate points to make a good plot. Usually it will use about 200 points, but you can adjust that to whatever number you want.

P:= plot(x^2, x= 0..3):
ExportMatrix(
     "C:/Users/Carl/desktop/Filename.txt", op([1,1], P),
     target= delimited, delimiter= " "
);

The code op([1,1], P) is what extracts the data matrix from the plot.

I am not sure what you mean by indenting the negative numbers. If the below is not it, let me know.

S:= map(x-> sprintf("%a", lhs(x)), v1):
S:= map(StringTools:-PadRight, S, max(length~(S))):
seq(printf("%s = %3.3f\n", S[k], rhs(v1[k])), k= 1..numelems(v1));

eta[p2]   = 0.260
eta[p3]   = 0.113
eta[p4]   = -0.013
eta[p5]   = 0.215
eta[p6]   = -0.189
eta[phi2] = 0.020
eta[phi3] = 0.063
eta[phi4] = -0.014
eta[phi5] = -0.414
eta[phi6] = 0.067
mu[p]     = 0.466
mu[phi]   = -0.169
tau[p3]   = 0.000
tau[p4]   = 0.000
w[1]      = 0.023
w[2]      = -0.447
w[3]      = -0.110
w[4]      = 0.035
w[5]      = 0.445