@k20057 5 

Before understanding the code you need to understand the combinatorial concepts of (proper) compositions and weak compositions. If you don't, see this Wikipedia page.

The command combinat:-composition(n+k, k) returns a set of all proper compositions of n+k of length k, each composition being represented by a list. Each composition of n+k needs to be converted to a weak composition of n by subtracting 1 from each element (because what you want is a weak composition of 8). The expression C-~1 subtracts 1 from every element of C. (The ~ makes a binary operator apply to every element of its operands.) The command seq(C-~1, C= L) iterates that over every element of L. The square brackets around the seq command turn the result into a list, which is then called Compositions. So, this is a list of all weak compositions of n of length k.

The expression nops(Compositions) returns the number of elements of Compositions. Apparently you already know what rand(a..b) does because you used it in your original Question.

The line ()-> Compositions[Rand()] is the return value of the procedure. This return value is itself a procedure created with the -> operator. The () to the left of -> indicates that this procedure takes no arguments (or has no parameters). The Rand() generates a random integer in the appropriate range. The expression L[n] selects the nth member of list L.

@k20057 5 

Here it is.

RandomCompositions:= proc(n::posint, k::posint)

(* Generates a procedure which selects uniformly at random a k-tuple of nonnegative
integers that sums to n. Note that this has a memory requirement of
k*binomial(n+k-1, k-1).                                                                                                            *)

local
     C,
     Compositions:= [seq(C-~1, C= combinat:-composition(n+k, k))],
     Rand:= rand(1..nops(Compositions))
;
     ()-> Compositions[Rand()]
end proc:

R:= RandomCompositions(8,6):
R();

R();

@liushunyi You wrote:

In my opinion, data is just the file name of input and ofile is the file name of output. Could we remove the declarations data := "read-data-compute-charpoly.dat" and ofile := "read-data-compute-charpoly.poly"?

Joe intended that you would change the assignments to data and ofile to be the names of your actual files.

The error is caused by extra white space characters at the end of your input file. It can be corrected like this. The line

if line=0 or line="" or line[1]="#" then

should be changed to

if line = 0 or line::string and StringTools:-TrimRight(line) = "" or line[1] = "#" then

I did this, and the program ran without error, producing the correct output file.

@Aysan 

Referring to the same system of five equations under discussion in this Answer, you can set the relerr and abserr options like this:

Sol1:= dsolve({Sys1, ICs1},
     numeric, maxfun= 0,
     method= rosenbrock,
     abserr= 1e-10, relerr= 1e-6,
     range= 0..10
):

You should also set, globally,

Digits:=  15:

increasing it from its default value of 10. That's the "sweet spot" for simultaneous accuracy and speed.

Please read the help page ?dsolve,rosenbrock .

@Aysan

According to ?dsolve,rosenbrock, the default values for both relative and absolute tolerance (options relerr and abserr, respectively) are both 10^(-6). Where did you see 10^(-3)?

@Aysan 

Consider the following piecewise expression of three pieces:

f:= piecewise(
     t < 1,  #operand 1
     t,      #operand 2
     t < 2,  #operand 3
     t^2,    #operand 4
     t >= 2, #operand 5
     t^3     #operand 6
);

So, the Boolean conditions tend to be the odd-numbered operands, and the algebraic "pieces" tend to be the even-numbered operands.

@Kitonum I like this better than my solution with diff.

Since dsolve returns multiple solutions, the odetest call needs to be

odetest({sol}, ode);

Can an expression that contains a RootOf where the RootOf variable appears as a limit of integration of an undoable integral really be considered to be a "solution" of the ODE?

@Aysan 

Since eq1 is independent of x, the answer is just a constant multiple of the original answer, the constant being int(sin(Pi*x), x= 0..1), which is 2/Pi.

@smith_alpha 

As shown by Kitonum, my procedure does not work when the base variable appears in the exponent.

@Aysan 

f:= t-> piecewise(t < 10, 1-t, t):
ode:= diff(y(t),t$2) - y(t)^2 = f(t):
Sol:= dsolve({ode, y(0)=0, D(y)(0)=0}, numeric):
plots:-odeplot(Sol, [t, y(t)], t= 0..12);

@Kitonum 

Your code fails for many special cases, such as when the variable appears to power 1, when the variable appears multiple times in the expression, and when the variable does not appear in the expression.

@liushunyi To obtain just the sequence do

seq(coeff(p, x, k), k= degree(p,x)..0, -1);

or

PolynomialTools:-CoefficientList(p, x, termorder= reverse)[];

In general, if L is a list or set, then L[] is the corresponding sequence. The same thing can be achieved with op(L), but the operation is so common that that just clutters up the code.

@Mac Dude Indeed, I'd recommend that the OP use the interface command that you gave because "I don't want to click on 'properties' each time I execute my worksheet."

Please give us the commands that you are using to do the export.