Can you give an example please?

If the dependent variable is not changing much over an interval of width h, that corresponds with its derivative being small. So, you can make the events like this (assuming that de is in the form with the derivative on the lhs):

events= [[abs(rhs(de))  - 1e-4, halt]]

@spradlig See ?list or ?set (it's the same page), especially the 7th and 8th paragraphs of "Description":

Appending an element x to a list L is done by [op(L), x].  Inserting an element x to a set S is done using the union operator S union {x}.
Replacing the i-th element of a list L by x can be done by subsop(i=x, L). Deleting the i-th element of a list L is subsop(i=NULL, L). Deleting an element x from a set S is done using the minus operator S minus {x}.  Note that inserting and deleting from a list or set causes a new list or set to be created rather than modifying the original in-place.

@abbeykabir Enclose the equations with curly braces { } when passing to fsolve. It probably will use Newton's Method and the Jacobian.

Here's an example. I can't cut-and-paste your equations, so I generated random similar examples.

'randpoly([y[k] $ k= 1..5], degree= 3, sparse) + randpoly([y[k] $ k= 1..5], degree= 1)' $ 5;

fsolve({%});

@abbeykabir Enclose the equations with curly braces { } when passing to fsolve. It probably will use Newton's Method and the Jacobian.

Here's an example. I can't cut-and-paste your equations, so I generated random similar examples.

'randpoly([y[k] $ k= 1..5], degree= 3, sparse) + randpoly([y[k] $ k= 1..5], degree= 1)' $ 5;

fsolve({%});

Your #2 is quite clever, but the syntax requires it to be 2..nops(L2) rather than 2:nops(L2).

Your #2 is quite clever, but the syntax requires it to be 2..nops(L2) rather than 2:nops(L2).

Indeed, this more general rule will work:

applyrule(a::algebraic^c::algebraic*b::algebraic^c::algebraic= (a*b)^c, %);

Yet combine(%, power, symbolic) does not work. (Why?)

Indeed, this more general rule will work:

applyrule(a::algebraic^c::algebraic*b::algebraic^c::algebraic= (a*b)^c, %);

Yet combine(%, power, symbolic) does not work. (Why?)

Yet there are the symbolic options to simplify, sqrt, and combine(..., power). One should be able to achieve what the Asker wants under those options, throwing caution to the wind, so to speak. I can't understand why the following does not work:

1/sqrt(x)*1/sqrt(y);
combine(%, power, symbolic);

But the applyrule seems to work easily.

Yet there are the symbolic options to simplify, sqrt, and combine(..., power). One should be able to achieve what the Asker wants under those options, throwing caution to the wind, so to speak. I can't understand why the following does not work:

1/sqrt(x)*1/sqrt(y);
combine(%, power, symbolic);

But the applyrule seems to work easily.

If we add a prettyprinting procedure to Markiyan's code, we can get output that looks more like the original formula without altering any of the computational aspects:

`print/binomial`:= (n,k)-> 'C'[n]^k:

theo(a1,a2,N1,N2);

If we add a prettyprinting procedure to Markiyan's code, we can get output that looks more like the original formula without altering any of the computational aspects:

`print/binomial`:= (n,k)-> 'C'[n]^k:

theo(a1,a2,N1,N2);

Since expressions are assumed to be equated to 0 by solve, and since the solution variables are assumed to be all the indeterminates, you could simply do

solve(map(coeffs, {f,g}, [x,y]));

Since expressions are assumed to be equated to 0 by solve, and since the solution variables are assumed to be all the indeterminates, you could simply do

solve(map(coeffs, {f,g}, [x,y]));