Advanced Search

mmcdara

7284 Reputation

18 Badges

8 years, 259 days
Contact mmcdara

MaplePrimes Activity


These are answers submitted by mmcdara

inconsistent unknowns...

mmcdara 7284
December 31 2024
3 5

Let P a multidimensional parameter.
You have 

phis := {phin = Fn(G(xi), P) , n=1..7 }

and you want to solve the system 

{Fn(G(xi), P) = 0 , n=1..7 }

wrt to P.

To get P you use the command

unknowns := indets(rhs~(phis), name)

which indeed returns all the names P contains, but also the name (xi) of the independent variable G depends upon.
So the error you get:
Error, (in solve) cannot solve expressions with diff(G(xi), xi) for xi

To get rid of it write

unknowns := indets(rhs~(phis), name) minus {xi};   # this is the "true" P

Now you can compute COEFFS, which is a sequence of 16 solutions:

COEFFS := solve(rhs~(phis), unknowns)

whattype(COEFFS);
                            exprseq

 numelems([COEFFS]);
                               16

Some of them depend on G(xi), which might not (?) be suitable.
If you want to get rid of them execute

GfreeCOEFFS := remove(has, [COEFFS], xi):
numelems(%);
                               9

Remark:

phis contains 7 relations while the parameter vector P has dimension 16. Which mean that 9 parameters are likely to get arbitrary values.

Maybe a better way would be to solve the system  { Fn(G(xi), Q) = 0 , n=1..7 } with Q being a subvector parameter of P with size  7.
Ideally you should chose the 7 parameters which are the most relevant to the use you have in mind.

Is I'm not mistaken......

mmcdara 7284
December 30 2024
2 3

and ifI understand corrrectly your problem, here is a solution (which is likely not optimal in terms of computation time nor memory use) .
It is fast for N=8 but I didn't checked it for much larger values of N.

Round_Robin_mmcdara.mw

Maple Worksheet - Error

Failed to load the worksheet /maplenet/convert/round_robin_mmcdara.mw .



 

A way...

mmcdara 7284
December 27 2024
4 7


key: use this rewritting rule

isolate(m+1/(diff(G(xi), xi))=Phi, diff(G(xi), xi));
                       d              1   
                      ---- G(xi) = -------
                       dxi         Phi - m

and collect wrt Phi :

> 

restart
> 

with(PDEtools):
> 

with(LinearAlgebra):
> 

> 

with(SolveTools):
> 

undeclare(prime)

`There is no more prime differentiation variable; all derivatives will be displayed as indexed functions`

 (1)
> 

PDEtools:-declare(Omega(x, t)); 1; PDEtools:-declare(U(xi)); 1; PDEtools:-declare(u(x, y, z, t)); 1; PDEtools:-declare(Q(xi))

Omega(x, t)*`will now be displayed as`*Omega

  

U(xi)*`will now be displayed as`*U

  

u(x, y, z, t)*`will now be displayed as`*u

  

Q(xi)*`will now be displayed as`*Q

 (2)
> 

tr := {t = tau, x = (-ZETA*c[3]-tau*c[4]-`Υ`*c[2]+xi)/c[1], y = `Υ`, z = ZETA, u(x, y, z, t) = U(xi)}:
> 

pde1 := diff(u(x, y, z, t), `$`(x, 3), z)-4*(diff(u(x, y, z, t), x, t))+4*(diff(u(x, y, z, t), x))*(diff(u(x, y, z, t), x, z))+2*(diff(u(x, y, z, t), `$`(x, 2)))*(diff(u(x, y, z, t), z))+3*(diff(u(x, y, z, t), `$`(y, 2))) = 0:
> 

NULL
> 

L1 := PDEtools:-dchange(tr, pde1, [xi, `Υ`, ZETA, tau, U]):
> 

map(int, L1, xi):
> 

ode := %:
> 

F := sum(a[i]*(m+1/(diff(G(xi), xi)))^i, i = -1 .. 1):
> 

D1 := diff(F, xi):
> 

S := diff(G(xi), `$`(xi, 2)) = -(2*m*mu+lambda)*(diff(G(xi), xi))-mu:
> 

E1 := subs(S, D1):
> 

D2 := diff(E1, xi):
> 

E2 := subs(S, D2):
> 

D3 := diff(E2, xi):
> 

E3 := subs(S, D3):
> 

``
> 

> 

K := U(xi) = F:
> 

K1 := diff(U(xi), xi) = E1:
> 

K2 := diff(U(xi), `$`(xi, 2)) = E2:
> 

K3 := diff(U(xi), `$`(xi, 3)) = E3:
> 

``
> 

L := eval(ode, {K, K1, K2, K3}):
> 

> 

# rewritting rule

RR := isolate(m+1/(diff(G(xi), xi))=Phi, diff(G(xi), xi));

diff(G(xi), xi) = 1/(Phi-m)

 (3)
> 

# Apply RR and collect wrt Phi

subs(RR, L):
normal(%):
PhiN := collect(numer(lhs(%)), Phi):
PhiD := denom(lhs(%%));

Phi^4

 (4)
> 



with(LargeExpressions):

LLE := collect(PhiN, Phi, Veil[phi] ):
LLE / PhiD = 0;

(3*Phi^8*phi[1]+6*Phi^7*phi[2]+Phi^6*phi[3]+Phi^5*phi[4]-Phi^4*phi[5]-Phi^3*phi[6]+Phi^2*phi[7]-6*Phi*phi[8]+3*phi[9])/Phi^4 = 0

 (5)
> 

# phi[i] coefficients

print~( [ seq( phi[i] = simplify(Unveil[phi](phi[i]), size), i=1..LastUsed[phi] ) ] ):

phi[1] = mu^2*a[1]*c[1]^2*c[3]*(2*mu*c[1]+a[1])

  

phi[2] = lambda*mu*a[1]*c[1]^2*c[3]*(2*mu*c[1]+a[1])

  

phi[3] = -8*a[1]*(mu*c[3]*(m^2*mu^2+m*mu*lambda-(7/8)*lambda^2)*c[1]^3+(3/4)*((m^2*a[1]+a[-1])*mu^2+m*mu*lambda*a[1]-(1/2)*lambda^2*a[1])*c[3]*c[1]^2+(1/2)*mu*c[1]*c[4]-(3/8)*mu*c[2]^2)

  

phi[4] = -8*lambda*a[1]*(c[3]*(m^2*mu^2+m*mu*lambda-(1/8)*lambda^2)*c[1]^3+(3/4)*(m*lambda*a[1]+mu*(m^2*a[1]+2*a[-1]))*c[3]*c[1]^2+(1/2)*c[1]*c[4]-(3/8)*c[2]^2)

  

phi[5] = -2*((m^2*a[1]+a[-1])*mu+m*lambda*a[1])*c[3]*(m^2*mu^2+m*mu*lambda-(1/2)*lambda^2)*c[1]^3-3*((m^4*a[1]^2+4*m^2*a[-1]*a[1]+a[-1]^2)*mu^2+2*m*lambda*a[1]*(m^2*a[1]+2*a[-1])*mu+lambda^2*a[1]*(m^2*a[1]-2*a[-1]))*c[3]*c[1]^2-4*c[4]*((m^2*a[1]+a[-1])*mu+m*lambda*a[1])*c[1]+3*((m^2*a[1]+a[-1])*mu+m*lambda*a[1])*c[2]^2

  

phi[6] = -8*lambda*(c[3]*(m^2*mu^2+m*mu*lambda-(1/8)*lambda^2)*c[1]^3+(3/2)*(m*lambda*a[1]+mu*(m^2*a[1]+(1/2)*a[-1]))*c[3]*c[1]^2+(1/2)*c[1]*c[4]-(3/8)*c[2]^2)*a[-1]

  

phi[7] = -8*(mu^2*c[1]^2*c[3]*(mu*c[1]+(3/4)*a[1])*m^4+2*lambda*(mu*c[1]+(3/4)*a[1])*mu*c[3]*c[1]^2*m^3+((1/8)*mu*c[3]*c[1]^3*lambda^2+(3/4)*c[3]*(lambda^2*a[1]+mu^2*a[-1])*c[1]^2+(1/2)*mu*c[1]*c[4]-(3/8)*mu*c[2]^2)*m^2+(3/4)*lambda*(-(7/6)*c[1]^3*c[3]*lambda^2+mu*a[-1]*c[1]^2*c[3]+(2/3)*c[1]*c[4]-(1/2)*c[2]^2)*m-(3/8)*lambda^2*a[-1]*c[1]^2*c[3])*a[-1]

  

phi[8] = lambda*m*c[1]^2*(2*m^2*mu*c[1]+2*lambda*m*c[1]+a[-1])*c[3]*(m*mu+lambda)*a[-1]

  

phi[9] = m^2*c[1]^2*(2*m^2*mu*c[1]+2*lambda*m*c[1]+a[-1])*c[3]*(m*mu+lambda)^2*a[-1]

 (6)
 

 

Download factoring_mmcdara.mw

All answer is problem and machine depend...

mmcdara 7284
December 27 2024
3 4


Nevertheless here is a real example based on the Grid package and a specific strategy (read carefully the text).
 

> 

restart
> 

with(Grid):


Run N_sim simulations on N_noeud nodes (can be larger than Numnodes() )

Here I want to solve an ode system to determine the temporal evolution of X(t ; P) for different
values of a parameter vector P.

Those values are gathered  in a matrix named Ver_plan_exp. of dimensions N_sim by N_P where
N_P is the size of P.

Ver_plan_exp is considered as made of N_blocs blocks of equal sizes T_blocs by N_P.

Those N_blocs blocs are run independently on N_noeud nodes: in case N_blocs is larger than N_noeud
the first N_noeud blocs are simulated and once all have ended their job another bunch of N_noeud blocks
is sent to them dor simulation. The process keeps going on until all the N_blocs blocks are simulated.

This very simple strategy (Grid:-Wait command) is acceptable as the simulation times for each bloc are roughly
the same (so I use here some kind of synchronous strategy).
When it's not the case it's better to use an asynchronous strategy: as soon as one of the active node ens its task its is
charged with a new vlock without waiting for the other to send and ending signal. This strategy is suitable when the
different blocs lead to significantly different execution times.
Finally, you can notice there is no master node here as all the N_noeud nodes do exactly the same task.
I could have use a special node, let's say node 0, to monitor the tasks of the remaining  N_noeud-1 
slave nodes. This master node then have the task to distribute the simulation and to gather ther results these latters
deliver.
At this point it's worth saying that, in my case, the procedure run on each node (SIMULER_MP) generates a huge
matrix made of  T_bloc rows and hundreds of columns. If each node has to return such a matrix to a possible
master node in order this latter aggregate them, this will have a cost in terms of communication time and memory
used. This is why I prefered not to use a master node here.

Of course the choices between  synchronous vs  asynchronous strategy and  all-slaves vs one-master strategy are
problem dependent.
But they also depend on the OS you use and of what the term node means: are the nodes physical or logical.
The strategy I present here is the simplest one, and it proved to be the most efficient for my problem on my machine
(a 24 nodes PC (12 physical) with 64 Gb memory under Windows 10... for information I got better performances using
Windows XP which didn't interfere with the Grid commands, see below).
At last, the value of N_bloc may impact the efficiency: too small block sizes mean more communications and thus
penalize the computational time, but this can limit the memory consumption. So here again a "good" choice will depend on
your machine and the OS you use.


The key procedure here is SIMULER_MP whose aim is to sole the ode system for all the parameter (row) vectors
P defined by the submatrix plan_bloc = Ver_plan_exp[bloc].

This procedure is declared this way
SIMULER_MP := proc(
                    num_bloc,
                    entree,
                    sortie,
                    Drag_sol,
                    tps_decollement,
                    Ver_piecewise_va,
                    mur,
                    tps_max,
                    ICI
              )

As I said above the solutions computed by this procedure are gathered in a big matrix named RESULT.
I found that the most efficient solution was, after the construction of this matrix, to save it in binary format
before ending procedure SIMULER_MP :

save RESULT, cat(ICI, "/tps_", num_bloc, ".m");
end proc:


In the code below you see that how the nine arguments of SIMULER_MP are passed in the Grid:-Run command.

Once all the simulations are done it may be interesting to kill the processes associated to all the nodes but node 0 (which
is the default number associated to your worksheet) to freeze the memory.
The way to do that depends on the OS you use.

It remains to read all the m files and to gather their content, if necessary,in a single m file.


Last advices:

(1)  it may be interesting to use a number of nodes lower than the number Grid:-NumNodes returns if you
       want to assume some light tasks aside (mailing, text edititong..., things like that),
(2)   keep an aye on the monitor panel to see the nodes are correctly loaded.

        (when the code below was ran under Window XP the task monitor displayed 20 nodes with a 100% charge,
        when ranunder Windows 10 it this was no longer the case because Windows 10 [and presumablyt higher versions]
        use its own managing process to load the nodes. The result was that, on the same machine, the computational time
        was 20-25% higher with Windows 10 than with Windows XP).


Good luck.

> 

Grid:-NumNodes();  # returns the number of nodes (starts from 0 by default)

N_noeud   := 20;
N_blocs   := 20;
N_sim     := 10^6
T_blocs   := N_sim / N_blocs;

Grid:-Setup(numnodes=N_noeud):
Grid:-Set  ('SIMULER_MP'):

printf("%s   départ\n", StringTools:-FormatTime("%H:%M:%S"));

for m from 1 to N_blocs/N_noeud do
  for noeud from 0 to N_noeud-1 do
    num_bloc  := (m-1)*N_noeud+noeud;
    bloc      := 1 + num_bloc*T_blocs..(num_bloc+1)*T_blocs;
    plan_bloc := Ver_plan_exp[bloc]:

    #-------------------------------------------------
    Grid:-Run(
      noeud,                  # target node
      SIMULER_MP,             # procedure to run
      [                        
        num_bloc,             #
        plan_bloc,            #
        [1],                  # list of arguments
        eval(Drag_sol),       #
        tps_decollement,      #
        [Ver_piecewise_va],   #
        mur,                  #
        tps_max               #
        ICI                   #
      ]                       
    ):
    #-------------------------------------------------
  end do;

  Grid:-Wait(seq(0..N_noeud-1));

  printf("%s   Fin groupe %2d\n", StringTools:-FormatTime("%H:%M:%S"), m);
end do;

printf("%s   Fin\n", StringTools:-FormatTime("%H:%M:%S"));


 

Download Example.mw

Bad use if a code written for something ...

mmcdara 7284
December 25 2024
1 3


The piece of code

	
W__points := plottools:-getdata(P)[1, -1]:
t_grid := convert(W__points[..,1], list):
x_grid := [seq(-2..2, 4/N)]:
> 	
T, X := map(mul, [selectremove(has, [op(expand(solnum))], t)])[]:
ST := unapply(eval(T, W(t)=w), w)~(W__points[.., 2]):
SX := evalf(unapply(X, x)~(x_grid)):
STX := Matrix(N$2, (it, ix) -> ST[it]*SX[ix])

I sent you days ago was written to handle T3 expressions of the form F(x)*G(W(t)), which is no longer the case here.

Here is an update for your current T3 expression (in your code you gave w two different meanings: one was a parameter, the other a surrogate name for W(t), which is quite confusing ; even if I no longer use this surrogate in the code below, I thought it would be better to avoid confusion in any future  expressions like eval(..., W(t)=w) and then arbitrarily renamed w as psi).

 

> 

restart;
> 

currentdir(kernelopts(':-homedir')):
> 

randomize():
> 

local gamma:
> 

T3 := (B[1]*(tanh(2*n^2*(delta^2-psi)*k*t/((k*n-1)*(k*n+1))+x)-1))^(1/(2*n))*exp(I*(-k*x+psi*t+delta*W(t)-delta^2*t));

(B[1]*(tanh(2*n^2*(delta^2-psi)*k*t/((k*n-1)*(k*n+1))+x)-1))^((1/2)/n)*exp(I*(-k*x+psi*t+delta*W(t)-delta^2*t))

 (1)
> 

params := {B[1]=1, n=2, delta=2, psi=1, k=3 }:

insert numerical values

> 

solnum :=subs(params, T3):
> 

N := 100:
use Finance in:
  Wiener := WienerProcess():
  P := PathPlot(Wiener(t), t = 0..10, timesteps = N, replications = 1):
end use:
> 

W__points := plottools:-getdata(P)[1, -1];
t_grid := convert(W__points[..,1], list):
x_grid := [seq(-2..2, 4/N)]:

`#msub(mi("W"),mi("points"))` := Vector(4, {(1) = ` 101 x 2 `*Matrix, (2) = `Data Type: `*float[8], (3) = `Storage: `*rectangular, (4) = `Order: `*Fortran_order})

 (2)
> 

T3_param := eval(T3, params):

STX := Matrix(
         (N+1)$2
         , (it, ix) ->
           evalf(
             eval(
               eval( T3_param, W(t)=W__points[it, 2])
               , {t=t_grid[it], x=x_grid[ix]}
             )
           )
         );

STX := Vector(4, {(1) = ` 101 x 101 `*Matrix, (2) = `Data Type: `*anything, (3) = `Storage: `*rectangular, (4) = `Order: `*Fortran_order})

 (3)
> 

plots:-matrixplot(Re~(STX));

 
> 

# Export STX as you use to do
> 

 


 

Download S5_mmcdara.mw

output=Array...

mmcdara 7284
December 24 2024
1 0

The simplest way, IMO, to get this table  is to use the option output=Array in dsolve/numeric:

when    := [seq(0..50, 1)]:
sol     := dsolve(`union`(sys, ic), numeric, output=Array(when)):
sol[1];  # to verify the order of the outputs
results := round~(sol[2][1]):

printf("%-10s %-15s %-15s %-15s\n", "Day", "Infected", "Recovered", "Susceptible");
printf("---------------------------------------------\n");

fmt := "%-10d %-15d %-15d %-15d\n":
map(i -> printf(fmt, entries(results[i+1], nolist)), when):

Another way...

mmcdara 7284
December 23 2024
4 1 
y := substring(x, 1..1)

To get the value of y:

x := asdf:
a := 3:
y := substring(x, 1..1);
                               a
eval(y)
                               3

For information only...

mmcdara 7284
December 22 2024
0 2

This doesn't answer your question, just to say that geom3d provides an elegant way (IMO) to determine if a point is on a line.

restart
with(geom3d):
line(L, [3+2*alpha, 1+6*alpha, 4-5*alpha], alpha):

point(A, 9, 19, -11):
IsOnObject(A, L)
                              true
point(B, 9, 1, -11):
IsOnObject(B, L)
                             false

Use similar triangles...

mmcdara 7284
December 22 2024
3 0


Plus an animation to illustrate the demonstration
Similar_Triangles_Animation.mw

A less brute force...

mmcdara 7284
December 20 2024
2 1

Here is a much faster (still quite brute) way to get the result

> 

restart:
> 

N := n -> add(n[i]*10^(i-1), i=1..6):
> 

P := [ seq(ListTools:-Rotate([seq(n[i], i=1..6)], j), j=0..5) ]

[[n[1], n[2], n[3], n[4], n[5], n[6]], [n[2], n[3], n[4], n[5], n[6], n[1]], [n[3], n[4], n[5], n[6], n[1], n[2]], [n[4], n[5], n[6], n[1], n[2], n[3]], [n[5], n[6], n[1], n[2], n[3], n[4]], [n[6], n[1], n[2], n[3], n[4], n[5]]]

 (1)
> 

tstart := time():

NP   := N~(P):
Z    := NP[2..-1]:
ZP   := combinat:-permute(Z):
sols := NULL:
for m from 1 to numelems(ZP) do
  eval(isolve({seq(NP[1]*(k+1) = ZP[m][k], k=1..5)}), _Z1=1);
  if max(rhs~(%)) < 10 then
    sols := sols, %
  end if;
end do:
sols;

time()-tstart;

{n[1] = 7, n[2] = 5, n[3] = 8, n[4] = 2, n[5] = 4, n[6] = 1}

  

.502

 (2)
> 

N([eval(seq(n[i], i=1..6), sols)]);
% *~ [$1..6];
 

142857

  

[142857, 285714, 428571, 571428, 714285, 857142]

 (3)
> 

# vv code

tstart := time():

P:=combinat:-permute([seq](0..9),6): # nops(%);
> 

for L in P do
> 

  n:=parse(cat(L[]));
> 

  ok:=true;
> 

  for k from 2 to 6 do
> 

    n1:=n*k;
> 

    L1:=convert(n1,base,10):
> 

    if {L[]} <> {L1[]} then ok:=false; break fi;
> 

  od;
> 

if ok then print(n, x23456=[seq(k*n,k=2..6)]) fi;
> 

od:

time()-tstart
> 

 

142857, x23456 = [285714, 428571, 571428, 714285, 857142]

  

20.158

 (4)
 

 

Download Improvement.mw

Here it is...

mmcdara 7284
December 20 2024
2 5

After a lot of work as you didn't provide all the necessary details:

Note there is likely an error in the paper you present: the transformation should probably be 

u = -2 * ln(fxx) instead of u = 2 * ln(fxx)

, as this screen capture seems to confirm
 

> 

restart

 

Emulation of the Hirota derivation operator

> 

alias(F=F(x, t), G=G(x, t))

F, G

 (1)
> 

with(PDEtools):
undeclare(prime):

`There is no more prime differentiation variable; all derivatives will be displayed as indexed functions`

 (2)
> 

ND := proc(F, G, U)
  local v, w, f, g, a:
  v := op(F):
  if v[1] in U then w := -v[1] else w := v[1] end if:
  if v[2] in U then w := w, -v[2] else w := w, v[2] end if:
  f := op(0, F):
  g := op(0, G):
  a := diff(f(w)*g(v), U);
  convert(subs([w]=~[v], a), diff)
end proc:


Verify if
         ND(F, F, [x, t]) + ND(F, F, [x$6])
is indeed equal to the lhs of equation (3.2) in the excerpt of the paper you present.

> 

ND(f, f, [x, t]) + ND(f, f, [x$6]);

2*f*(diff(diff(f, t), x))-2*(diff(f, x))*(diff(f, t))+2*f*(diff(diff(diff(diff(diff(diff(f, x), x), x), x), x), x))-12*(diff(f, x))*(diff(diff(diff(diff(diff(f, x), x), x), x), x))+30*(diff(diff(f, x), x))*(diff(diff(diff(diff(f, x), x), x), x))-20*(diff(diff(diff(f, x), x), x))^2

 (3)

 

Sawada-Kotera equation

> 

alias(u=u(x, t), f=f(x, t))

F, G, u, w, f

 (4)
> 

SK__u := diff(u, t) + 45*u^2*diff(u, x) - 15*diff(u, x)*diff(u, x$2) - 15*u*diff(u, x$3) + diff(u, x$5) = 0

diff(u, t)+45*u^2*(diff(u, x))-15*(diff(u, x))*(diff(diff(u, x), x))-15*u*(diff(diff(diff(u, x), x), x))+diff(diff(diff(diff(diff(u, x), x), x), x), x) = 0

 (5)
> 

SK__f := eval(SK__u, u=diff(f, x$2))

diff(diff(diff(f, t), x), x)+45*(diff(diff(f, x), x))^2*(diff(diff(diff(f, x), x), x))-15*(diff(diff(diff(f, x), x), x))*(diff(diff(diff(diff(f, x), x), x), x))-15*(diff(diff(f, x), x))*(diff(diff(diff(diff(diff(f, x), x), x), x), x))+diff(diff(diff(diff(diff(diff(diff(f, x), x), x), x), x), x), x) = 0

 (6)
> 

ISK__f := map(Int, lhs(SK__f), x) = 0

Int(diff(diff(diff(f, t), x), x), x)+Int(45*(diff(diff(f, x), x))^2*(diff(diff(diff(f, x), x), x)), x)+Int(-15*(diff(diff(diff(f, x), x), x))*(diff(diff(diff(diff(f, x), x), x), x)), x)+Int(-15*(diff(diff(f, x), x))*(diff(diff(diff(diff(diff(f, x), x), x), x), x)), x)+Int(diff(diff(diff(diff(diff(diff(diff(f, x), x), x), x), x), x), x), x) = 0

 (7)
> 

ISK__fv := value(ISK__f);

diff(diff(f, t), x)+15*(diff(diff(f, x), x))^3-(15/2)*(diff(diff(diff(f, x), x), x))^2+int(-15*(diff(diff(f, x), x))*(diff(diff(diff(diff(diff(f, x), x), x), x), x)), x)+diff(diff(diff(diff(diff(diff(f, x), x), x), x), x), x) = 0

 (8)


The problem here is the uneval integral term.
Let Z its integrand:

> 

Z := IntegrationTools:-GetIntegrand( select(has, [op(lhs(ISK__fv))], int)[] );

-15*(diff(diff(f, x), x))*(diff(diff(diff(diff(diff(f, x), x), x), x), x))

 (9)


Use two by part integrations to compute its integral:

> 

use IntegrationTools in
  ``(Int(Z, x)) = Parts(Int(Z, x), diff(f, x$2));
  lhs(%) = expand(op(1, rhs(%)) + Parts(Expand(op(2, rhs(%))), diff(f, x$3)));
  % + %%;
  map(Expand, %);
  Reductor := expand~(% /~ 2);
end use;

``(Int(-15*(diff(diff(f, x), x))*(diff(diff(diff(diff(diff(f, x), x), x), x), x)), x)) = -15*(diff(diff(f, x), x))*(diff(diff(diff(diff(f, x), x), x), x))-(Int(-15*(diff(diff(diff(f, x), x), x))*(diff(diff(diff(diff(f, x), x), x), x)), x))

  

``(Int(-15*(diff(diff(f, x), x))*(diff(diff(diff(diff(diff(f, x), x), x), x), x)), x)) = -15*(diff(diff(f, x), x))*(diff(diff(diff(diff(f, x), x), x), x))+15*(diff(diff(diff(f, x), x), x))^2-15*(Int((diff(diff(diff(f, x), x), x))*(diff(diff(diff(diff(f, x), x), x), x)), x))

  

2*``(Int(-15*(diff(diff(f, x), x))*(diff(diff(diff(diff(diff(f, x), x), x), x), x)), x)) = -30*(diff(diff(f, x), x))*(diff(diff(diff(diff(f, x), x), x), x))+15*(diff(diff(diff(f, x), x), x))^2-15*(Int((diff(diff(diff(f, x), x), x))*(diff(diff(diff(diff(f, x), x), x), x)), x))-(Int(-15*(diff(diff(diff(f, x), x), x))*(diff(diff(diff(diff(f, x), x), x), x)), x))

  

2*``(-15*(Int((diff(diff(f, x), x))*(diff(diff(diff(diff(diff(f, x), x), x), x), x)), x))) = -30*(diff(diff(f, x), x))*(diff(diff(diff(diff(f, x), x), x), x))+15*(diff(diff(diff(f, x), x), x))^2

  

-15*(Int((diff(diff(f, x), x))*(diff(diff(diff(diff(diff(f, x), x), x), x), x)), x)) = -15*(diff(diff(f, x), x))*(diff(diff(diff(diff(f, x), x), x), x))+(15/2)*(diff(diff(diff(f, x), x), x))^2

 (10)


Use this relation to reduce ISK__f  and  get an expression free from integrals:

> 

eval(IntegrationTools:-Expand(ISK__f), Reductor);

ISK__fv := value(%);

Int(diff(diff(diff(f, t), x), x), x)+45*(Int((diff(diff(f, x), x))^2*(diff(diff(diff(f, x), x), x)), x))-15*(Int((diff(diff(diff(f, x), x), x))*(diff(diff(diff(diff(f, x), x), x), x)), x))-15*(diff(diff(f, x), x))*(diff(diff(diff(diff(f, x), x), x), x))+(15/2)*(diff(diff(diff(f, x), x), x))^2+Int(diff(diff(diff(diff(diff(diff(diff(f, x), x), x), x), x), x), x), x) = 0

  

diff(diff(f, t), x)+15*(diff(diff(f, x), x))^3-15*(diff(diff(f, x), x))*(diff(diff(diff(diff(f, x), x), x), x))+diff(diff(diff(diff(diff(diff(f, x), x), x), x), x), x) = 0

 (11)


Use the transformation F = -2*ln(f)

> 

alias(F = F(x, t)):

ISK__F := normal(eval(ISK__fv, f=alpha*ln(F))):

ISK__F := eval(numer(lhs(ISK__F)), alpha=-2) = 0;

-2*(diff(diff(diff(diff(diff(diff(F, x), x), x), x), x), x))*F^5+12*(diff(diff(diff(diff(diff(F, x), x), x), x), x))*(diff(F, x))*F^4-30*(diff(diff(diff(diff(F, x), x), x), x))*(diff(diff(F, x), x))*F^4+20*(diff(diff(diff(F, x), x), x))^2*F^4-2*(diff(diff(F, t), x))*F^5+2*(diff(F, t))*(diff(F, x))*F^4 = 0

 (12)
> 

# As the highest derivation degree wrt x is 6, the D-operator will contain
# ND(F, F, [x$6]):

x6 := `#msubsup(mo("D"),mo("x"),mo("6"))` = ND(F, F, [x$6]);

ToRewrite := diff(F, x$6):
RewriteAs := isolate(x6, ToRewrite);
Rewritten := simplify(subs( eval(RewriteAs, alpha=2), lhs(ISK__F)));

`#msubsup(mo("D"),mo("x"),mo("6"))` = 2*F*(diff(diff(diff(diff(diff(diff(F, x), x), x), x), x), x))-12*(diff(diff(diff(diff(diff(F, x), x), x), x), x))*(diff(F, x))+30*(diff(diff(diff(diff(F, x), x), x), x))*(diff(diff(F, x), x))-20*(diff(diff(diff(F, x), x), x))^2

  

diff(diff(diff(diff(diff(diff(F, x), x), x), x), x), x) = -(1/2)*(30*(diff(diff(diff(diff(F, x), x), x), x))*(diff(diff(F, x), x))-20*(diff(diff(diff(F, x), x), x))^2-12*(diff(diff(diff(diff(diff(F, x), x), x), x), x))*(diff(F, x))-`#msubsup(mo("D"),mo("x"),mo("6"))`)/F

  

-F^4*(2*(diff(diff(F, t), x))*F-2*(diff(F, t))*(diff(F, x))+`#msubsup(mo("D"),mo("x"),mo("6"))`)

 (13)
> 

# Rewrite the terms containing the second derivative of F wrt x and t

xt := `#msubsup(mo("D"),mrow(mo("x"),mo("t")))` = ND(F, F, [x, t]);

ToRewrite := diff(F, [x, t]):
RewriteAs := isolate(xt, ToRewrite);
Rewritten := collect(simplify(expand(subs(RewriteAs, Rewritten))), [F, `#msubsup(mo("D"),mo("x"),mo("2"))`])

`#msubsup(mo("D"),mrow(mo("x"),mo("t")))` = 2*(diff(diff(F, t), x))*F-2*(diff(F, t))*(diff(F, x))

  

diff(diff(F, t), x) = -(1/2)*(-2*(diff(F, t))*(diff(F, x))-`#msubsup(mo("D"),mrow(mo("x"),mo("t")))`)/F

  

(-`#msubsup(mo("D"),mo("x"),mo("6"))`-`#msubsup(mo("D"),mrow(mo("x"),mo("t")))`)*F^4

 (14)

 

Hirota*form*of*SK*equation*(alpha = -2)

 

> 

SK_H := ``(eval(-Rewritten, F=1))(F.F)=0

(``(`#msubsup(mo("D"),mrow(mo("x"),mo("t")))`+`#msubsup(mo("D"),mo("x"),mo("6"))`))(F.F) = 0

 (15)

 

"To get a relation closer to the las one in the excerpt of the paper you present we need to prove this; D[xt] = D[x] @ D[t] = D[t] @ D[x]"

 

> 

`#msub(mo("D"),mrow(mo("x"),mo("t")))` = ND(F, G, [x, t]);


# some algebraic transformations
     relation_x := `#msub(mo("D"),mo("x"))` = ND(F, G, [x]):
     alias(f=f(x, t), g=g(x, t)):
     equivalences := {diff(F, x)=f, diff(G, x)=g}:

     eval(relation_x, equivalences):
     map(z -> sign(z)*ND(op(sign(z)*z), [t]), [op(rhs(%))]):
     add(%):

`#msub(mo("D"),mo("x"))`@`#msub(mo("D"),mo("t"))` = eval(%, (rhs=lhs)~(equivalences));


# some algebraic transformations
     relation_t := `#msub(mo("D"),mo("t"))` = ND(F, G, [t]):
     equivalences := {diff(F, t)=f, diff(G, t)=g}:

     eval(relation_t, equivalences):
     map(z -> sign(z)*ND(op(sign(z)*z), [x]), [op(rhs(%))]):
     add(%):

`#msub(mo("D"),mo("t"))`@`#msub(mo("D"),mo("x"))` = eval(%, (rhs=lhs)~(equivalences))

`#msub(mo("D"),mrow(mo("x"),mo("t")))` = (diff(diff(F, t), x))*G-(diff(F, x))*(diff(G, t))-(diff(F, t))*(diff(G, x))+F*(diff(diff(G, t), x))

  

`@`(`#msub(mo("D"),mo("x"))`, `#msub(mo("D"),mo("t"))`) = (diff(diff(F, t), x))*G-(diff(F, x))*(diff(G, t))-(diff(F, t))*(diff(G, x))+F*(diff(diff(G, t), x))

  

`@`(`#msub(mo("D"),mo("t"))`, `#msub(mo("D"),mo("x"))`) = (diff(diff(F, t), x))*G-(diff(F, x))*(diff(G, t))-(diff(F, t))*(diff(G, x))+F*(diff(diff(G, t), x))

 (16)


Then

> 

SK_H := eval((15), { lhs(xt) = `#msub(mo("D"),mo("x"))` * `#msub(mo("D"),mo("t"))`, F=f })

(``(`#msub(mo("D"),mo("t"))`*`#msub(mo("D"),mo("x"))`+`#msubsup(mo("D"),mo("x"),mo("6"))`))(f.f) = 0

 (17)
 

 

Download Sawada_Kotera_(Hirota_form).mw

'known' option...

mmcdara 7284
December 19 2024
3 3

@rcorless already answered you.

Nevertheless, if for some personal reason you prefer to solve a first system with unknown functions x(t), y(t) and z(t), and then plug the solutions within a second system with unknown functions P(t), Q(t) and R(t), you can use the known option of dsolve/numeric (look to the help page).
This is a very powerful feeature which is worth being known as it can be use in several other situations.

Here is how to use it

> 

restart:
> 

eq1 := diff(x(t), t)-(1/6)*(6*x(t)^3*y(t)+(2*y(t)^2-2)*x(t)^2+3*y(t)*(z(t)-2)*x(t)-2*y(t)^2+2)*sqrt(3) = 0;
> 

                         
> 

eq2 := diff(y(t), t)-(1/6)*(y(t)-1)*sqrt(3)*(y(t)+1)*(6*x(t)^2+2*y(t)*x(t)+3*z(t)-2) = 0;
> 

                                    
> 

eq3 := diff(z(t), t)-(1/3)*z(t)*sqrt(3)*(6*y(t)*x(t)^2+2*x(t)*y(t)^2+3*z(t)*y(t)-2*x(t)-3*y(t)) = 0;
> 

 

> 

sys := eq1, eq2, eq3:
> 

ics := x(0) = -0.01, y(0) = .99, z(0) = 0.01;
> 

 

> 

sol := dsolve({ics, sys}, type = numeric);

diff(x(t), t)-(1/6)*(6*x(t)^3*y(t)+(2*y(t)^2-2)*x(t)^2+3*y(t)*(z(t)-2)*x(t)-2*y(t)^2+2)*3^(1/2) = 0

  

diff(y(t), t)-(1/6)*(y(t)-1)*3^(1/2)*(y(t)+1)*(6*x(t)^2+2*y(t)*x(t)+3*z(t)-2) = 0

  

diff(z(t), t)-(1/3)*z(t)*3^(1/2)*(6*y(t)*x(t)^2+2*x(t)*y(t)^2+3*z(t)*y(t)-2*x(t)-3*y(t)) = 0

  

x(0) = -0.1e-1, y(0) = .99, z(0) = 0.1e-1

  

proc (x_rkf45) local _res, _dat, _vars, _solnproc, _xout, _ndsol, _pars, _n, _i; option `Copyright (c) 2000 by Waterloo Maple Inc. All rights reserved.`; if 1 < nargs then error "invalid input: too many arguments" end if; _EnvDSNumericSaveDigits := Digits; Digits := 15; if _EnvInFsolve = true then _xout := evalf[_EnvDSNumericSaveDigits](x_rkf45) else _xout := evalf(x_rkf45) end if; _dat := Array(1..4, {(1) = proc (_xin) local _xout, _dtbl, _dat, _vmap, _x0, _y0, _val, _dig, _n, _ne, _nd, _nv, _pars, _ini, _par, _i, _j, _k, _src; option `Copyright (c) 2002 by Waterloo Maple Inc. All rights reserved.`; table( [( "complex" ) = false ] ) _xout := _xin; _pars := []; _dtbl := array( 1 .. 4, [( 1 ) = (array( 1 .. 24, [( 1 ) = (datatype = float[8], order = C_order, storage = rectangular), ( 2 ) = (datatype = float[8], order = C_order, storage = rectangular), ( 3 ) = ([0, 0, 0, Array(1..0, 1..2, {}, datatype = float[8], order = C_order)]), ( 4 ) = (Array(1..54, {(1) = 3, (2) = 3, (3) = 0, (4) = 0, (5) = 0, (6) = 0, (7) = 1, (8) = 0, (9) = 0, (10) = 0, (11) = 0, (12) = 0, (13) = 0, (14) = 0, (15) = 0, (16) = 0, (17) = 0, (18) = 1, (19) = 30000, (20) = 0, (21) = 0, (22) = 1, (23) = 4, (24) = 0, (25) = 1, (26) = 15, (27) = 1, (28) = 0, (29) = 1, (30) = 3, (31) = 3, (32) = 0, (33) = 1, (34) = 0, (35) = 0, (36) = 0, (37) = 0, (38) = 0, (39) = 0, (40) = 0, (41) = 0, (42) = 0, (43) = 1, (44) = 0, (45) = 0, (46) = 0, (47) = 0, (48) = 0, (49) = 0, (50) = 50, (51) = 1, (52) = 0, (53) = 0, (54) = 0}, datatype = integer[8])), ( 5 ) = (Array(1..28, {(1) = .0, (2) = 0.10e-5, (3) = .0, (4) = 0.500001e-14, (5) = .0, (6) = .19449453019010343, (7) = .0, (8) = 0.10e-5, (9) = .0, (10) = .0, (11) = .0, (12) = .0, (13) = 1.0, (14) = .0, (15) = .49999999999999, (16) = .0, (17) = 1.0, (18) = 1.0, (19) = .0, (20) = .0, (21) = 1.0, (22) = 1.0, (23) = .0, (24) = .0, (25) = 0.10e-14, (26) = .0, (27) = .0, (28) = .0}, datatype = float[8], order = C_order)), ( 6 ) = (Array(1..3, {(1) = -0.1e-1, (2) = .99, (3) = 0.1e-1}, datatype = float[8], order = C_order)), ( 7 ) = ([Array(1..4, 1..7, {(1, 1) = .0, (1, 2) = .203125, (1, 3) = .3046875, (1, 4) = .75, (1, 5) = .8125, (1, 6) = .40625, (1, 7) = .8125, (2, 1) = 0.6378173828125e-1, (2, 2) = .0, (2, 3) = .279296875, (2, 4) = .27237892150878906, (2, 5) = -0.9686851501464844e-1, (2, 6) = 0.1956939697265625e-1, (2, 7) = .5381584167480469, (3, 1) = 0.31890869140625e-1, (3, 2) = .0, (3, 3) = -.34375, (3, 4) = -.335235595703125, (3, 5) = .2296142578125, (3, 6) = .41748046875, (3, 7) = 11.480712890625, (4, 1) = 0.9710520505905151e-1, (4, 2) = .0, (4, 3) = .40350341796875, (4, 4) = 0.20297467708587646e-1, (4, 5) = -0.6054282188415527e-2, (4, 6) = -0.4770040512084961e-1, (4, 7) = .77858567237854}, datatype = float[8], order = C_order), Array(1..6, 1..6, {(1, 1) = .0, (1, 2) = .0, (1, 3) = .0, (1, 4) = .0, (1, 5) = .0, (1, 6) = 1.0, (2, 1) = .25, (2, 2) = .0, (2, 3) = .0, (2, 4) = .0, (2, 5) = .0, (2, 6) = 1.0, (3, 1) = .1875, (3, 2) = .5625, (3, 3) = .0, (3, 4) = .0, (3, 5) = .0, (3, 6) = 2.0, (4, 1) = .23583984375, (4, 2) = -.87890625, (4, 3) = .890625, (4, 4) = .0, (4, 5) = .0, (4, 6) = .2681884765625, (5, 1) = .1272735595703125, (5, 2) = -.5009765625, (5, 3) = .44921875, (5, 4) = -0.128936767578125e-1, (5, 5) = .0, (5, 6) = 0.626220703125e-1, (6, 1) = -0.927734375e-1, (6, 2) = .626220703125, (6, 3) = -.4326171875, (6, 4) = .1418304443359375, (6, 5) = -0.861053466796875e-1, (6, 6) = .3131103515625}, datatype = float[8], order = C_order), Array(1..6, {(1) = .0, (2) = .386, (3) = .21, (4) = .63, (5) = 1.0, (6) = 1.0}, datatype = float[8], order = C_order), Array(1..6, {(1) = .25, (2) = -.1043, (3) = .1035, (4) = -0.362e-1, (5) = .0, (6) = .0}, datatype = float[8], order = C_order), Array(1..6, 1..5, {(1, 1) = .0, (1, 2) = .0, (1, 3) = .0, (1, 4) = .0, (1, 5) = .0, (2, 1) = 1.544, (2, 2) = .0, (2, 3) = .0, (2, 4) = .0, (2, 5) = .0, (3, 1) = .9466785280815533, (3, 2) = .25570116989825814, (3, 3) = .0, (3, 4) = .0, (3, 5) = .0, (4, 1) = 3.3148251870684886, (4, 2) = 2.896124015972123, (4, 3) = .9986419139977808, (4, 4) = .0, (4, 5) = .0, (5, 1) = 1.2212245092262748, (5, 2) = 6.019134481287752, (5, 3) = 12.537083329320874, (5, 4) = -.687886036105895, (5, 5) = .0, (6, 1) = 1.2212245092262748, (6, 2) = 6.019134481287752, (6, 3) = 12.537083329320874, (6, 4) = -.687886036105895, (6, 5) = 1.0}, datatype = float[8], order = C_order), Array(1..6, 1..5, {(1, 1) = .0, (1, 2) = .0, (1, 3) = .0, (1, 4) = .0, (1, 5) = .0, (2, 1) = -5.6688, (2, 2) = .0, (2, 3) = .0, (2, 4) = .0, (2, 5) = .0, (3, 1) = -2.4300933568337584, (3, 2) = -.20635991570891224, (3, 3) = .0, (3, 4) = .0, (3, 5) = .0, (4, 1) = -.10735290581452621, (4, 2) = -9.594562251021896, (4, 3) = -20.470286148096154, (4, 4) = .0, (4, 5) = .0, (5, 1) = 7.496443313968615, (5, 2) = -10.246804314641219, (5, 3) = -33.99990352819906, (5, 4) = 11.708908932061595, (5, 5) = .0, (6, 1) = 8.083246795922411, (6, 2) = -7.981132988062785, (6, 3) = -31.52159432874373, (6, 4) = 16.319305431231363, (6, 5) = -6.0588182388340535}, datatype = float[8], order = C_order), Array(1..3, 1..5, {(1, 1) = .0, (1, 2) = .0, (1, 3) = .0, (1, 4) = .0, (1, 5) = .0, (2, 1) = 10.126235083446911, (2, 2) = -7.487995877607633, (2, 3) = -34.800918615557414, (2, 4) = -7.9927717075687275, (2, 5) = 1.0251377232956207, (3, 1) = -.6762803392806898, (3, 2) = 6.087714651678606, (3, 3) = 16.43084320892463, (3, 4) = 24.767225114183653, (3, 5) = -6.5943891257167815}, datatype = float[8], order = C_order)]), ( 9 ) = ([Array(1..3, {(1) = .1, (2) = .1, (3) = .1}, datatype = float[8], order = C_order), Array(1..3, {(1) = .0, (2) = .0, (3) = .0}, datatype = float[8], order = C_order), Array(1..3, {(1) = .0, (2) = .0, (3) = .0}, datatype = float[8], order = C_order), Array(1..3, {(1) = .0, (2) = .0, (3) = .0}, datatype = float[8], order = C_order), Array(1..3, {(1) = .0, (2) = .0, (3) = .0}, datatype = float[8], order = C_order), Array(1..3, 1..3, {(1, 1) = .0, (1, 2) = .0, (1, 3) = .0, (2, 1) = .0, (2, 2) = .0, (2, 3) = .0, (3, 1) = .0, (3, 2) = .0, (3, 3) = .0}, datatype = float[8], order = C_order), Array(1..3, 1..3, {(1, 1) = .0, (1, 2) = .0, (1, 3) = .0, (2, 1) = .0, (2, 2) = .0, (2, 3) = .0, (3, 1) = .0, (3, 2) = .0, (3, 3) = .0}, datatype = float[8], order = C_order), Array(1..3, 1..6, {(1, 1) = .0, (1, 2) = .0, (1, 3) = .0, (1, 4) = .0, (1, 5) = .0, (1, 6) = .0, (2, 1) = .0, (2, 2) = .0, (2, 3) = .0, (2, 4) = .0, (2, 5) = .0, (2, 6) = .0, (3, 1) = .0, (3, 2) = .0, (3, 3) = .0, (3, 4) = .0, (3, 5) = .0, (3, 6) = .0}, datatype = float[8], order = C_order), Array(1..3, {(1) = 0, (2) = 0, (3) = 0}, datatype = integer[8]), Array(1..3, {(1) = .0, (2) = .0, (3) = .0}, datatype = float[8], order = C_order), Array(1..3, {(1) = .0, (2) = .0, (3) = .0}, datatype = float[8], order = C_order), Array(1..3, {(1) = .0, (2) = .0, (3) = .0}, datatype = float[8], order = C_order), Array(1..3, {(1) = .0, (2) = .0, (3) = .0}, datatype = float[8], order = C_order), Array(1..3, {(1) = .0, (2) = .0, (3) = .0}, datatype = float[8], order = C_order)]), ( 8 ) = ([Array(1..3, {(1) = -0.1e-1, (2) = .99, (3) = 0.1e-1}, datatype = float[8], order = C_order), Array(1..3, {(1) = .0, (2) = .0, (3) = .0}, datatype = float[8], order = C_order), Array(1..3, {(1) = 0.28547973179495756e-1, (2) = 0.1142722829694645e-1, (3) = -0.1697010265031224e-1}, datatype = float[8], order = C_order), 0, 0]), ( 11 ) = (Array(1..6, 0..3, {(1, 1) = .0, (1, 2) = .0, (1, 3) = .0, (2, 0) = .0, (2, 1) = .0, (2, 2) = .0, (2, 3) = .0, (3, 0) = .0, (3, 1) = .0, (3, 2) = .0, (3, 3) = .0, (4, 0) = .0, (4, 1) = .0, (4, 2) = .0, (4, 3) = .0, (5, 0) = .0, (5, 1) = .0, (5, 2) = .0, (5, 3) = .0, (6, 0) = .0, (6, 1) = .0, (6, 2) = .0, (6, 3) = .0}, datatype = float[8], order = C_order)), ( 10 ) = ([proc (N, X, Y, YP) option `[Y[1] = x(t), Y[2] = y(t), Y[3] = z(t)]`; YP[1] := 1.73205080756888*Y[1]^3*Y[2]+.288675134594813*(2*Y[2]^2-2)*Y[1]^2+.866025403784440*Y[2]*(Y[3]-2)*Y[1]-.577350269189627*Y[2]^2+.577350269189627; YP[2] := .288675134594813*(Y[2]-1)*(Y[2]+1)*(6*Y[1]^2+2*Y[1]*Y[2]+3*Y[3]-2); YP[3] := .577350269189627*Y[3]*(6*Y[1]^2*Y[2]+2*Y[1]*Y[2]^2+3*Y[2]*Y[3]-2*Y[1]-3*Y[2]); 0 end proc, -1, 0, 0, 0, 0, 0, 0]), ( 13 ) = (), ( 12 ) = (), ( 15 ) = ("rkf45"), ( 14 ) = ([0, 0]), ( 18 ) = ([]), ( 19 ) = (0), ( 16 ) = ([0, 0, 0, []]), ( 17 ) = ([proc (N, X, Y, YP) option `[Y[1] = x(t), Y[2] = y(t), Y[3] = z(t)]`; YP[1] := 1.73205080756888*Y[1]^3*Y[2]+.288675134594813*(2*Y[2]^2-2)*Y[1]^2+.866025403784440*Y[2]*(Y[3]-2)*Y[1]-.577350269189627*Y[2]^2+.577350269189627; YP[2] := .288675134594813*(Y[2]-1)*(Y[2]+1)*(6*Y[1]^2+2*Y[1]*Y[2]+3*Y[3]-2); YP[3] := .577350269189627*Y[3]*(6*Y[1]^2*Y[2]+2*Y[1]*Y[2]^2+3*Y[2]*Y[3]-2*Y[1]-3*Y[2]); 0 end proc, -1, 0, 0, 0, 0, 0, 0]), ( 22 ) = (0), ( 23 ) = (0), ( 20 ) = ([]), ( 21 ) = (0), ( 24 ) = (0) ] )) ] ); _y0 := Array(0..3, {(1) = 0., (2) = -0.1e-1, (3) = .99}); _vmap := array( 1 .. 3, [( 1 ) = (1), ( 2 ) = (2), ( 3 ) = (3) ] ); _x0 := _dtbl[1][5][5]; _n := _dtbl[1][4][1]; _ne := _dtbl[1][4][3]; _nd := _dtbl[1][4][4]; _nv := _dtbl[1][4][16]; if not type(_xout, 'numeric') then if member(_xout, ["start", "left", "right"]) then if _Env_smart_dsolve_numeric = true or _dtbl[1][4][10] = 1 then if _xout = "left" then if type(_dtbl[2], 'table') then return _dtbl[2][5][1] end if elif _xout = "right" then if type(_dtbl[3], 'table') then return _dtbl[3][5][1] end if end if end if; return _dtbl[1][5][5] elif _xout = "method" then return _dtbl[1][15] elif _xout = "storage" then return evalb(_dtbl[1][4][10] = 1) elif _xout = "leftdata" then if not type(_dtbl[2], 'array') then return NULL else return eval(_dtbl[2]) end if elif _xout = "rightdata" then if not type(_dtbl[3], 'array') then return NULL else return eval(_dtbl[3]) end if elif _xout = "enginedata" then return eval(_dtbl[1]) elif _xout = "enginereset" then _dtbl[2] := evaln(_dtbl[2]); _dtbl[3] := evaln(_dtbl[3]); return NULL elif _xout = "initial" then return procname(_y0[0]) elif _xout = "laxtol" then return _dtbl[`if`(member(_dtbl[4], {2, 3}), _dtbl[4], 1)][5][18] elif _xout = "numfun" then return `if`(member(_dtbl[4], {2, 3}), _dtbl[_dtbl[4]][4][18], 0) elif _xout = "parameters" then return [seq(_y0[_n+_i], _i = 1 .. nops(_pars))] elif _xout = "initial_and_parameters" then return procname(_y0[0]), [seq(_y0[_n+_i], _i = 1 .. nops(_pars))] elif _xout = "last" then if _dtbl[4] <> 2 and _dtbl[4] <> 3 or _x0-_dtbl[_dtbl[4]][5][1] = 0. then error "no information is available on last computed point" else _xout := _dtbl[_dtbl[4]][5][1] end if elif _xout = "function" then if _dtbl[1][4][33]-2. = 0 then return eval(_dtbl[1][10], 1) else return eval(_dtbl[1][10][1], 1) end if elif _xout = "map" then return copy(_vmap) elif type(_xin, `=`) and type(rhs(_xin), 'list') and member(lhs(_xin), {"initial", "parameters", "initial_and_parameters"}) then _ini, _par := [], []; if lhs(_xin) = "initial" then _ini := rhs(_xin) elif lhs(_xin) = "parameters" then _par := rhs(_xin) elif select(type, rhs(_xin), `=`) <> [] then _par, _ini := selectremove(type, rhs(_xin), `=`) elif nops(rhs(_xin)) < nops(_pars)+1 then error "insufficient data for specification of initial and parameters" else _par := rhs(_xin)[-nops(_pars) .. -1]; _ini := rhs(_xin)[1 .. -nops(_pars)-1] end if; _xout := lhs(_xout); if _par <> [] then `dsolve/numeric/process_parameters`(_n, _pars, _par, _y0) end if; if _ini <> [] then `dsolve/numeric/process_initial`(_n-_ne, _ini, _y0, _pars, _vmap) end if; `dsolve/numeric/SC/reinitialize`(_dtbl, _y0, _n, procname, _pars); if _Env_smart_dsolve_numeric = true and type(_y0[0], 'numeric') and _dtbl[1][4][10] <> 1 then procname("right") := _y0[0]; procname("left") := _y0[0] end if; if _xout = "initial" then return [_y0[0], seq(_y0[_vmap[_i]], _i = 1 .. _n-_ne)] elif _xout = "parameters" then return [seq(_y0[_n+_i], _i = 1 .. nops(_pars))] else return [_y0[0], seq(_y0[_vmap[_i]], _i = 1 .. _n-_ne)], [seq(_y0[_n+_i], _i = 1 .. nops(_pars))] end if elif _xin = "eventstop" then if _nv = 0 then error "this solution has no events" end if; _i := _dtbl[4]; if _i <> 2 and _i <> 3 then return 0 end if; if _dtbl[_i][4][10] = 1 and assigned(_dtbl[5-_i]) and _dtbl[_i][4][9] < 100 and 100 <= _dtbl[5-_i][4][9] then _i := 5-_i; _dtbl[4] := _i; _j := round(_dtbl[_i][4][17]); return round(_dtbl[_i][3][1][_j, 1]) elif 100 <= _dtbl[_i][4][9] then _j := round(_dtbl[_i][4][17]); return round(_dtbl[_i][3][1][_j, 1]) else return 0 end if elif _xin = "eventstatus" then if _nv = 0 then error "this solution has no events" end if; _i := [selectremove(proc (a) options operator, arrow; _dtbl[1][3][1][a, 7] = 1 end proc, {seq(_j, _j = 1 .. round(_dtbl[1][3][1][_nv+1, 1]))})]; return ':-enabled' = _i[1], ':-disabled' = _i[2] elif _xin = "eventclear" then if _nv = 0 then error "this solution has no events" end if; _i := _dtbl[4]; if _i <> 2 and _i <> 3 then error "no events to clear" end if; if _dtbl[_i][4][10] = 1 and assigned(_dtbl[5-_i]) and _dtbl[_i][4][9] < 100 and 100 < _dtbl[5-_i][4][9] then _dtbl[4] := 5-_i; _i := 5-_i end if; if _dtbl[_i][4][9] < 100 then error "no events to clear" elif _nv < _dtbl[_i][4][9]-100 then error "event error condition cannot be cleared" else _j := _dtbl[_i][4][9]-100; if irem(round(_dtbl[_i][3][1][_j, 4]), 2) = 1 then error "retriggerable events cannot be cleared" end if; _j := round(_dtbl[_i][3][1][_j, 1]); for _k to _nv do if _dtbl[_i][3][1][_k, 1] = _j then if _dtbl[_i][3][1][_k, 2] = 3 then error "range events cannot be cleared" end if; _dtbl[_i][3][1][_k, 8] := _dtbl[_i][3][1][_nv+1, 8] end if end do; _dtbl[_i][4][17] := 0; _dtbl[_i][4][9] := 0; if _dtbl[1][4][10] = 1 then if _i = 2 then try procname(procname("left")) catch: end try else try procname(procname("right")) catch: end try end if end if end if; return elif type(_xin, `=`) and member(lhs(_xin), {"eventdisable", "eventenable"}) then if _nv = 0 then error "this solution has no events" end if; if type(rhs(_xin), {('list')('posint'), ('set')('posint')}) then _i := {op(rhs(_xin))} elif type(rhs(_xin), 'posint') then _i := {rhs(_xin)} else error "event identifiers must be integers in the range 1..%1", round(_dtbl[1][3][1][_nv+1, 1]) end if; if select(proc (a) options operator, arrow; _nv < a end proc, _i) <> {} then error "event identifiers must be integers in the range 1..%1", round(_dtbl[1][3][1][_nv+1, 1]) end if; _k := {}; for _j to _nv do if member(round(_dtbl[1][3][1][_j, 1]), _i) then _k := `union`(_k, {_j}) end if end do; _i := _k; if lhs(_xin) = "eventdisable" then _dtbl[4] := 0; _j := [evalb(assigned(_dtbl[2]) and member(_dtbl[2][4][17], _i)), evalb(assigned(_dtbl[3]) and member(_dtbl[3][4][17], _i))]; for _k in _i do _dtbl[1][3][1][_k, 7] := 0; if assigned(_dtbl[2]) then _dtbl[2][3][1][_k, 7] := 0 end if; if assigned(_dtbl[3]) then _dtbl[3][3][1][_k, 7] := 0 end if end do; if _j[1] then for _k to _nv+1 do if _k <= _nv and not type(_dtbl[2][3][4][_k, 1], 'undefined') then userinfo(3, {'events', 'eventreset'}, `reinit #2, event code `, _k, ` to defined init `, _dtbl[2][3][4][_k, 1]); _dtbl[2][3][1][_k, 8] := _dtbl[2][3][4][_k, 1] elif _dtbl[2][3][1][_k, 2] = 0 and irem(iquo(round(_dtbl[2][3][1][_k, 4]), 32), 2) = 1 then userinfo(3, {'events', 'eventreset'}, `reinit #2, event code `, _k, ` to rate hysteresis init `, _dtbl[2][5][24]); _dtbl[2][3][1][_k, 8] := _dtbl[2][5][24] elif _dtbl[2][3][1][_k, 2] = 0 and irem(iquo(round(_dtbl[2][3][1][_k, 4]), 2), 2) = 0 then userinfo(3, {'events', 'eventreset'}, `reinit #2, event code `, _k, ` to initial init `, _x0); _dtbl[2][3][1][_k, 8] := _x0 else userinfo(3, {'events', 'eventreset'}, `reinit #2, event code `, _k, ` to fireinitial init `, _x0-1); _dtbl[2][3][1][_k, 8] := _x0-1 end if end do; _dtbl[2][4][17] := 0; _dtbl[2][4][9] := 0; if _dtbl[1][4][10] = 1 then procname(procname("left")) end if end if; if _j[2] then for _k to _nv+1 do if _k <= _nv and not type(_dtbl[3][3][4][_k, 2], 'undefined') then userinfo(3, {'events', 'eventreset'}, `reinit #3, event code `, _k, ` to defined init `, _dtbl[3][3][4][_k, 2]); _dtbl[3][3][1][_k, 8] := _dtbl[3][3][4][_k, 2] elif _dtbl[3][3][1][_k, 2] = 0 and irem(iquo(round(_dtbl[3][3][1][_k, 4]), 32), 2) = 1 then userinfo(3, {'events', 'eventreset'}, `reinit #3, event code `, _k, ` to rate hysteresis init `, _dtbl[3][5][24]); _dtbl[3][3][1][_k, 8] := _dtbl[3][5][24] elif _dtbl[3][3][1][_k, 2] = 0 and irem(iquo(round(_dtbl[3][3][1][_k, 4]), 2), 2) = 0 then userinfo(3, {'events', 'eventreset'}, `reinit #3, event code `, _k, ` to initial init `, _x0); _dtbl[3][3][1][_k, 8] := _x0 else userinfo(3, {'events', 'eventreset'}, `reinit #3, event code `, _k, ` to fireinitial init `, _x0+1); _dtbl[3][3][1][_k, 8] := _x0+1 end if end do; _dtbl[3][4][17] := 0; _dtbl[3][4][9] := 0; if _dtbl[1][4][10] = 1 then procname(procname("right")) end if end if else for _k in _i do _dtbl[1][3][1][_k, 7] := 1 end do; _dtbl[2] := evaln(_dtbl[2]); _dtbl[3] := evaln(_dtbl[3]); _dtbl[4] := 0; if _dtbl[1][4][10] = 1 then if _x0 <= procname("right") then try procname(procname("right")) catch: end try end if; if procname("left") <= _x0 then try procname(procname("left")) catch: end try end if end if end if; return elif type(_xin, `=`) and lhs(_xin) = "eventfired" then if not type(rhs(_xin), 'list') then error "'eventfired' must be specified as a list" end if; if _nv = 0 then error "this solution has no events" end if; if _dtbl[4] <> 2 and _dtbl[4] <> 3 then error "'direction' must be set prior to calling/setting 'eventfired'" end if; _i := _dtbl[4]; _val := NULL; if not assigned(_EnvEventRetriggerWarned) then _EnvEventRetriggerWarned := false end if; for _k in rhs(_xin) do if type(_k, 'integer') then _src := _k elif type(_k, 'integer' = 'anything') and type(evalf(rhs(_k)), 'numeric') then _k := lhs(_k) = evalf[max(Digits, 18)](rhs(_k)); _src := lhs(_k) else error "'eventfired' entry is not valid: %1", _k end if; if _src < 1 or round(_dtbl[1][3][1][_nv+1, 1]) < _src then error "event identifiers must be integers in the range 1..%1", round(_dtbl[1][3][1][_nv+1, 1]) end if; _src := {seq(`if`(_dtbl[1][3][1][_j, 1]-_src = 0., _j, NULL), _j = 1 .. _nv)}; if nops(_src) <> 1 then error "'eventfired' can only be set/queried for root-finding events and time/interval events" end if; _src := _src[1]; if _dtbl[1][3][1][_src, 2] <> 0. and _dtbl[1][3][1][_src, 2]-2. <> 0. then error "'eventfired' can only be set/queried for root-finding events and time/interval events" elif irem(round(_dtbl[1][3][1][_src, 4]), 2) = 1 then if _EnvEventRetriggerWarned = false then WARNING(`'eventfired' has no effect on events that retrigger`) end if; _EnvEventRetriggerWarned := true end if; if _dtbl[_i][3][1][_src, 2] = 0 and irem(iquo(round(_dtbl[_i][3][1][_src, 4]), 32), 2) = 1 then _val := _val, undefined elif type(_dtbl[_i][3][4][_src, _i-1], 'undefined') or _i = 2 and _dtbl[2][3][1][_src, 8] < _dtbl[2][3][4][_src, 1] or _i = 3 and _dtbl[3][3][4][_src, 2] < _dtbl[3][3][1][_src, 8] then _val := _val, _dtbl[_i][3][1][_src, 8] else _val := _val, _dtbl[_i][3][4][_src, _i-1] end if; if type(_k, `=`) then if _dtbl[_i][3][1][_src, 2] = 0 and irem(iquo(round(_dtbl[_i][3][1][_src, 4]), 32), 2) = 1 then error "cannot set event code for a rate hysteresis event" end if; userinfo(3, {'events', 'eventreset'}, `manual set event code `, _src, ` to value `, rhs(_k)); _dtbl[_i][3][1][_src, 8] := rhs(_k); _dtbl[_i][3][4][_src, _i-1] := rhs(_k) end if end do; return [_val] elif type(_xin, `=`) and lhs(_xin) = "direction" then if not member(rhs(_xin), {-1, 1, ':-left', ':-right'}) then error "'direction' must be specified as either '1' or 'right' (positive) or '-1' or 'left' (negative)" end if; _src := `if`(_dtbl[4] = 2, -1, `if`(_dtbl[4] = 3, 1, undefined)); _i := `if`(member(rhs(_xin), {1, ':-right'}), 3, 2); _dtbl[4] := _i; _dtbl[_i] := `dsolve/numeric/SC/IVPdcopy`(_dtbl[1], `if`(assigned(_dtbl[_i]), _dtbl[_i], NULL)); if 0 < _nv then for _j to _nv+1 do if _j <= _nv and not type(_dtbl[_i][3][4][_j, _i-1], 'undefined') then userinfo(3, {'events', 'eventreset'}, `reinit #4, event code `, _j, ` to defined init `, _dtbl[_i][3][4][_j, _i-1]); _dtbl[_i][3][1][_j, 8] := _dtbl[_i][3][4][_j, _i-1] elif _dtbl[_i][3][1][_j, 2] = 0 and irem(iquo(round(_dtbl[_i][3][1][_j, 4]), 32), 2) = 1 then userinfo(3, {'events', 'eventreset'}, `reinit #4, event code `, _j, ` to rate hysteresis init `, _dtbl[_i][5][24]); _dtbl[_i][3][1][_j, 8] := _dtbl[_i][5][24] elif _dtbl[_i][3][1][_j, 2] = 0 and irem(iquo(round(_dtbl[_i][3][1][_j, 4]), 2), 2) = 0 then userinfo(3, {'events', 'eventreset'}, `reinit #4, event code `, _j, ` to initial init `, _x0); _dtbl[_i][3][1][_j, 8] := _x0 else userinfo(3, {'events', 'eventreset'}, `reinit #4, event code `, _j, ` to fireinitial init `, _x0-2*_i+5.0); _dtbl[_i][3][1][_j, 8] := _x0-2*_i+5.0 end if end do end if; return _src elif _xin = "eventcount" then if _dtbl[1][3][1] = 0 or _dtbl[4] <> 2 and _dtbl[4] <> 3 then return 0 else return round(_dtbl[_dtbl[4]][3][1][_nv+1, 12]) end if else return "procname" end if end if; if _xout = _x0 then return [_x0, seq(evalf(_dtbl[1][6][_vmap[_i]]), _i = 1 .. _n-_ne)] end if; _i := `if`(_x0 <= _xout, 3, 2); if _xin = "last" and 0 < _dtbl[_i][4][9] and _dtbl[_i][4][9] < 100 then _dat := eval(_dtbl[_i], 2); _j := _dat[4][20]; return [_dat[11][_j, 0], seq(_dat[11][_j, _vmap[_i]], _i = 1 .. _n-_ne-_nd), seq(_dat[8][1][_vmap[_i]], _i = _n-_ne-_nd+1 .. _n-_ne)] end if; if not type(_dtbl[_i], 'array') then _dtbl[_i] := `dsolve/numeric/SC/IVPdcopy`(_dtbl[1], `if`(assigned(_dtbl[_i]), _dtbl[_i], NULL)); if 0 < _nv then for _j to _nv+1 do if _j <= _nv and not type(_dtbl[_i][3][4][_j, _i-1], 'undefined') then userinfo(3, {'events', 'eventreset'}, `reinit #5, event code `, _j, ` to defined init `, _dtbl[_i][3][4][_j, _i-1]); _dtbl[_i][3][1][_j, 8] := _dtbl[_i][3][4][_j, _i-1] elif _dtbl[_i][3][1][_j, 2] = 0 and irem(iquo(round(_dtbl[_i][3][1][_j, 4]), 32), 2) = 1 then userinfo(3, {'events', 'eventreset'}, `reinit #5, event code `, _j, ` to rate hysteresis init `, _dtbl[_i][5][24]); _dtbl[_i][3][1][_j, 8] := _dtbl[_i][5][24] elif _dtbl[_i][3][1][_j, 2] = 0 and irem(iquo(round(_dtbl[_i][3][1][_j, 4]), 2), 2) = 0 then userinfo(3, {'events', 'eventreset'}, `reinit #5, event code `, _j, ` to initial init `, _x0); _dtbl[_i][3][1][_j, 8] := _x0 else userinfo(3, {'events', 'eventreset'}, `reinit #5, event code `, _j, ` to fireinitial init `, _x0-2*_i+5.0); _dtbl[_i][3][1][_j, 8] := _x0-2*_i+5.0 end if end do end if end if; if _xin <> "last" then if 0 < 0 then if `dsolve/numeric/checkglobals`(op(_dtbl[1][14]), _pars, _n, _y0) then `dsolve/numeric/SC/reinitialize`(_dtbl, _y0, _n, procname, _pars, _i) end if end if; if _dtbl[1][4][7] = 0 then error "parameters must be initialized before solution can be computed" end if end if; _dat := eval(_dtbl[_i], 2); _dtbl[4] := _i; try _src := `dsolve/numeric/SC/IVPrun`(_dat, _xout) catch: userinfo(2, `dsolve/debug`, print(`Exception in solnproc:`, [lastexception][2 .. -1])); error end try; if _src = 0 and 100 < _dat[4][9] then _val := _dat[3][1][_nv+1, 8] else _val := _dat[11][_dat[4][20], 0] end if; if _src <> 0 or _dat[4][9] <= 0 then _dtbl[1][5][1] := _xout else _dtbl[1][5][1] := _val end if; if _i = 3 and _val < _xout then Rounding := -infinity; if _dat[4][9] = 1 then error "cannot evaluate the solution further right of %1, probably a singularity", evalf[8](_val) elif _dat[4][9] = 2 then error "cannot evaluate the solution further right of %1, maxfun limit exceeded (see ?dsolve,maxfun for details)", evalf[8](_val) elif _dat[4][9] = 3 then if _dat[4][25] = 3 then error "cannot evaluate the solution past the initial point, problem may be initially singular or improperly set up" else error "cannot evaluate the solution past the initial point, problem may be complex, initially singular or improperly set up" end if elif _dat[4][9] = 4 then error "cannot evaluate the solution further right of %1, accuracy goal cannot be achieved with specified 'minstep'", evalf[8](_val) elif _dat[4][9] = 5 then error "cannot evaluate the solution further right of %1, too many step failures, tolerances may be too loose for problem", evalf[8](_val) elif _dat[4][9] = 6 then error "cannot evaluate the solution further right of %1, cannot downgrade delay storage for problems with delay derivative order > 1, try increasing delaypts", evalf[8](_val) elif _dat[4][9] = 10 then error "cannot evaluate the solution further right of %1, interrupt requested", evalf[8](_val) elif 100 < _dat[4][9] then if _dat[4][9]-100 = _nv+1 then error "constraint projection failure on event at t=%1", evalf[8](_val) elif _dat[4][9]-100 = _nv+2 then error "index-1 and derivative evaluation failure on event at t=%1", evalf[8](_val) elif _dat[4][9]-100 = _nv+3 then error "maximum number of event iterations reached (%1) at t=%2", round(_dat[3][1][_nv+1, 3]), evalf[8](_val) else if _Env_dsolve_nowarnstop <> true then `dsolve/numeric/warning`(StringTools:-FormatMessage("cannot evaluate the solution further right of %1, event #%2 triggered a halt", evalf[8](_val), round(_dat[3][1][_dat[4][9]-100, 1]))) end if; Rounding := 'nearest'; _xout := _val end if else error "cannot evaluate the solution further right of %1", evalf[8](_val) end if elif _i = 2 and _xout < _val then Rounding := infinity; if _dat[4][9] = 1 then error "cannot evaluate the solution further left of %1, probably a singularity", evalf[8](_val) elif _dat[4][9] = 2 then error "cannot evaluate the solution further left of %1, maxfun limit exceeded (see ?dsolve,maxfun for details)", evalf[8](_val) elif _dat[4][9] = 3 then if _dat[4][25] = 3 then error "cannot evaluate the solution past the initial point, problem may be initially singular or improperly set up" else error "cannot evaluate the solution past the initial point, problem may be complex, initially singular or improperly set up" end if elif _dat[4][9] = 4 then error "cannot evaluate the solution further left of %1, accuracy goal cannot be achieved with specified 'minstep'", evalf[8](_val) elif _dat[4][9] = 5 then error "cannot evaluate the solution further left of %1, too many step failures, tolerances may be too loose for problem", evalf[8](_val) elif _dat[4][9] = 6 then error "cannot evaluate the solution further left of %1, cannot downgrade delay storage for problems with delay derivative order > 1, try increasing delaypts", evalf[8](_val) elif _dat[4][9] = 10 then error "cannot evaluate the solution further right of %1, interrupt requested", evalf[8](_val) elif 100 < _dat[4][9] then if _dat[4][9]-100 = _nv+1 then error "constraint projection failure on event at t=%1", evalf[8](_val) elif _dat[4][9]-100 = _nv+2 then error "index-1 and derivative evaluation failure on event at t=%1", evalf[8](_val) elif _dat[4][9]-100 = _nv+3 then error "maximum number of event iterations reached (%1) at t=%2", round(_dat[3][1][_nv+1, 3]), evalf[8](_val) else if _Env_dsolve_nowarnstop <> true then `dsolve/numeric/warning`(StringTools:-FormatMessage("cannot evaluate the solution further left of %1, event #%2 triggered a halt", evalf[8](_val), round(_dat[3][1][_dat[4][9]-100, 1]))) end if; Rounding := 'nearest'; _xout := _val end if else error "cannot evaluate the solution further left of %1", evalf[8](_val) end if end if; if _EnvInFsolve = true then _dig := _dat[4][26]; _dat[4][26] := _EnvDSNumericSaveDigits; _Env_dsolve_SC_native := true; if _dat[4][25] = 1 then _i := 1; _dat[4][25] := 2 else _i := _dat[4][25] end if; _val := `dsolve/numeric/SC/IVPval`(_dat, _xout, _src); _dat[4][25] := _i; _dat[4][26] := _dig; [_xout, seq(_val[_vmap[_i]], _i = 1 .. _n-_ne)] else Digits := _dat[4][26]; _val := `dsolve/numeric/SC/IVPval`(eval(_dat, 2), _xout, _src); [_xout, seq(_val[_vmap[_i]], _i = 1 .. _n-_ne)] end if end proc, (2) = Array(0..0, {}), (3) = [t, x(t), y(t), z(t)], (4) = []}); _vars := _dat[3]; _pars := map(rhs, _dat[4]); _n := nops(_vars)-1; _solnproc := _dat[1]; if not type(_xout, 'numeric') then if member(x_rkf45, ["start", 'start', "method", 'method', "left", 'left', "right", 'right', "leftdata", "rightdata", "enginedata", "eventstop", 'eventstop', "eventclear", 'eventclear', "eventstatus", 'eventstatus', "eventcount", 'eventcount', "laxtol", 'laxtol', "numfun", 'numfun', NULL]) then _res := _solnproc(convert(x_rkf45, 'string')); if 1 < nops([_res]) then return _res elif type(_res, 'array') then return eval(_res, 1) elif _res <> "procname" then return _res end if elif member(x_rkf45, ["last", 'last', "initial", 'initial', "parameters", 'parameters', "initial_and_parameters", 'initial_and_parameters', NULL]) then _xout := convert(x_rkf45, 'string'); _res := _solnproc(_xout); if _xout = "parameters" then return [seq(_pars[_i] = _res[_i], _i = 1 .. nops(_pars))] elif _xout = "initial_and_parameters" then return [seq(_vars[_i+1] = [_res][1][_i+1], _i = 0 .. _n), seq(_pars[_i] = [_res][2][_i], _i = 1 .. nops(_pars))] else return [seq(_vars[_i+1] = _res[_i+1], _i = 0 .. _n)] end if elif type(_xout, `=`) and member(lhs(_xout), ["initial", 'initial', "parameters", 'parameters', "initial_and_parameters", 'initial_and_parameters', NULL]) then _xout := convert(lhs(x_rkf45), 'string') = rhs(x_rkf45); if type(rhs(_xout), 'list') then _res := _solnproc(_xout) else error "initial and/or parameter values must be specified in a list" end if; if lhs(_xout) = "initial" then return [seq(_vars[_i+1] = _res[_i+1], _i = 0 .. _n)] elif lhs(_xout) = "parameters" then return [seq(_pars[_i] = _res[_i], _i = 1 .. nops(_pars))] else return [seq(_vars[_i+1] = [_res][1][_i+1], _i = 0 .. _n), seq(_pars[_i] = [_res][2][_i], _i = 1 .. nops(_pars))] end if elif type(_xout, `=`) and member(lhs(_xout), ["eventdisable", 'eventdisable', "eventenable", 'eventenable', "eventfired", 'eventfired', "direction", 'direction', NULL]) then return _solnproc(convert(lhs(x_rkf45), 'string') = rhs(x_rkf45)) elif _xout = "solnprocedure" then return eval(_solnproc) elif _xout = "sysvars" then return _vars end if; if procname <> unknown then return ('procname')(x_rkf45) else _ndsol; _ndsol := pointto(_dat[2][0]); return ('_ndsol')(x_rkf45) end if end if; try _res := _solnproc(_xout); [seq(_vars[_i+1] = _res[_i+1], _i = 0 .. _n)] catch: error end try end proc

 (1)
> 

eq4 := diff(R(t), t)-P(t)*Z(t)-(-2*(-Y(t)^2+2)*X(t)/sqrt(3)+sqrt(3)*(-2*X(t)^2+Z(t)+4/3)*Y(t))*R(t) = 0;
> 

eq5 := diff(Q(t), t)-(2/3)*R(t)+2*((1/3)*Y(t)+X(t))*P(t)/sqrt(3)-(-2*(-Y(t)^2+2)*X(t)/sqrt(3)+2*sqrt(3)*(X(t)^2-(1/2)*Z(t)-2/3)*X(t))*Q(t) = 0;
> 

eq6 := diff(P(t), t)+(1/2)*R(t)+2*sqrt(3)*X(t)*Q(t)+(2*(-Y(t)^2+2)*X(t)/sqrt(3)+sqrt(3)*(-2*X(t)^2+Z(t)+1)*Y(t))*P(t) = 0;
> 

 

> 

SYS := eq4, eq5, eq6:

diff(R(t), t)-P(t)*Z(t)-(-(2/3)*(-Y(t)^2+2)*X(t)*3^(1/2)+3^(1/2)*(-2*X(t)^2+Z(t)+4/3)*Y(t))*R(t) = 0

  

diff(Q(t), t)-(2/3)*R(t)+(2/3)*((1/3)*Y(t)+X(t))*P(t)*3^(1/2)-(-(2/3)*(-Y(t)^2+2)*X(t)*3^(1/2)+2*3^(1/2)*(X(t)^2-(1/2)*Z(t)-2/3)*X(t))*Q(t) = 0

  

diff(P(t), t)+(1/2)*R(t)+2*3^(1/2)*X(t)*Q(t)+((2/3)*(-Y(t)^2+2)*X(t)*3^(1/2)+3^(1/2)*(-2*X(t)^2+Z(t)+1)*Y(t))*P(t) = 0

 (2)
> 

indets({SYS});

{t, P(t), Q(t), R(t), X(t), Y(t), Z(t), diff(P(t), t), diff(Q(t), t), diff(R(t), t)}

 (3)
> 

# You didn't specicify ics ==> I use arbitrary ones

ICS := P(0) = rand(0. .. 1.)(), Q(0) = rand(0. .. 1.)(), R(0) = rand(0. .. 1.)()

P(0) = .2342493224, Q(0) = .1799302829, R(0) = .5137385362

 (4)
> 

X := proc(s)
  if s::numeric then
    eval(x(t), sol(s))
  else
    'procname'(s)
 end if:
end proc:



Y := proc(s)
  if s::numeric then
    eval(y(t), sol(s))
  else
    'procname'(s)
 end if:
end proc:



Z := proc(s)
  if s::numeric then
    eval(z(t), sol(s))
  else
    'procname'(s)
 end if:
end proc:


 
> 

SOL := dsolve({SYS, ICS}, type=numeric, known=[X(t), Y(t), Z(t)])

proc (x_rkf45) local _res, _dat, _vars, _solnproc, _xout, _ndsol, _pars, _n, _i; option `Copyright (c) 2000 by Waterloo Maple Inc. All rights reserved.`; if 1 < nargs then error "invalid input: too many arguments" end if; _EnvDSNumericSaveDigits := Digits; Digits := 15; if _EnvInFsolve = true then _xout := evalf[_EnvDSNumericSaveDigits](x_rkf45) else _xout := evalf(x_rkf45) end if; _dat := Array(1..4, {(1) = proc (_xin) local _xout, _dtbl, _dat, _vmap, _x0, _y0, _val, _dig, _n, _ne, _nd, _nv, _pars, _ini, _par, _i, _j, _k, _src; option `Copyright (c) 2002 by Waterloo Maple Inc. All rights reserved.`; table( [( "complex" ) = false ] ) _xout := _xin; _pars := []; _dtbl := array( 1 .. 4, [( 1 ) = (array( 1 .. 24, [( 1 ) = (datatype = float[8], order = C_order, storage = rectangular), ( 2 ) = (datatype = float[8], order = C_order, storage = rectangular), ( 3 ) = ([0, 0, 0, Array(1..0, 1..2, {}, datatype = float[8], order = C_order)]), ( 4 ) = (Array(1..54, {(1) = 3, (2) = 3, (3) = 0, (4) = 0, (5) = 0, (6) = 0, (7) = 1, (8) = 0, (9) = 0, (10) = 0, (11) = 0, (12) = 0, (13) = 0, (14) = 0, (15) = 0, (16) = 0, (17) = 0, (18) = 1, (19) = 30000, (20) = 0, (21) = 0, (22) = 1, (23) = 4, (24) = 0, (25) = 1, (26) = 15, (27) = 1, (28) = 0, (29) = 1, (30) = 3, (31) = 3, (32) = 0, (33) = 1, (34) = 0, (35) = 0, (36) = 0, (37) = 0, (38) = 0, (39) = 0, (40) = 0, (41) = 0, (42) = 0, (43) = 1, (44) = 0, (45) = 0, (46) = 0, (47) = 0, (48) = 0, (49) = 0, (50) = 50, (51) = 1, (52) = 0, (53) = 0, (54) = 0}, datatype = integer[8])), ( 5 ) = (Array(1..28, {(1) = .0, (2) = 0.10e-5, (3) = .0, (4) = 0.500001e-14, (5) = .0, (6) = 0.2581800197474211e-1, (7) = .0, (8) = 0.10e-5, (9) = .0, (10) = .0, (11) = .0, (12) = .0, (13) = 1.0, (14) = .0, (15) = .49999999999999, (16) = .0, (17) = 1.0, (18) = 1.0, (19) = .0, (20) = .0, (21) = 1.0, (22) = 1.0, (23) = .0, (24) = .0, (25) = 0.10e-14, (26) = .0, (27) = .0, (28) = .0}, datatype = float[8], order = C_order)), ( 6 ) = (Array(1..3, {(1) = .2342493224, (2) = .1799302829, (3) = .5137385362}, datatype = float[8], order = C_order)), ( 7 ) = ([Array(1..4, 1..7, {(1, 1) = .0, (1, 2) = .203125, (1, 3) = .3046875, (1, 4) = .75, (1, 5) = .8125, (1, 6) = .40625, (1, 7) = .8125, (2, 1) = 0.6378173828125e-1, (2, 2) = .0, (2, 3) = .279296875, (2, 4) = .27237892150878906, (2, 5) = -0.9686851501464844e-1, (2, 6) = 0.1956939697265625e-1, (2, 7) = .5381584167480469, (3, 1) = 0.31890869140625e-1, (3, 2) = .0, (3, 3) = -.34375, (3, 4) = -.335235595703125, (3, 5) = .2296142578125, (3, 6) = .41748046875, (3, 7) = 11.480712890625, (4, 1) = 0.9710520505905151e-1, (4, 2) = .0, (4, 3) = .40350341796875, (4, 4) = 0.20297467708587646e-1, (4, 5) = -0.6054282188415527e-2, (4, 6) = -0.4770040512084961e-1, (4, 7) = .77858567237854}, datatype = float[8], order = C_order), Array(1..6, 1..6, {(1, 1) = .0, (1, 2) = .0, (1, 3) = .0, (1, 4) = .0, (1, 5) = .0, (1, 6) = 1.0, (2, 1) = .25, (2, 2) = .0, (2, 3) = .0, (2, 4) = .0, (2, 5) = .0, (2, 6) = 1.0, (3, 1) = .1875, (3, 2) = .5625, (3, 3) = .0, (3, 4) = .0, (3, 5) = .0, (3, 6) = 2.0, (4, 1) = .23583984375, (4, 2) = -.87890625, (4, 3) = .890625, (4, 4) = .0, (4, 5) = .0, (4, 6) = .2681884765625, (5, 1) = .1272735595703125, (5, 2) = -.5009765625, (5, 3) = .44921875, (5, 4) = -0.128936767578125e-1, (5, 5) = .0, (5, 6) = 0.626220703125e-1, (6, 1) = -0.927734375e-1, (6, 2) = .626220703125, (6, 3) = -.4326171875, (6, 4) = .1418304443359375, (6, 5) = -0.861053466796875e-1, (6, 6) = .3131103515625}, datatype = float[8], order = C_order), Array(1..6, {(1) = .0, (2) = .386, (3) = .21, (4) = .63, (5) = 1.0, (6) = 1.0}, datatype = float[8], order = C_order), Array(1..6, {(1) = .25, (2) = -.1043, (3) = .1035, (4) = -0.362e-1, (5) = .0, (6) = .0}, datatype = float[8], order = C_order), Array(1..6, 1..5, {(1, 1) = .0, (1, 2) = .0, (1, 3) = .0, (1, 4) = .0, (1, 5) = .0, (2, 1) = 1.544, (2, 2) = .0, (2, 3) = .0, (2, 4) = .0, (2, 5) = .0, (3, 1) = .9466785280815533, (3, 2) = .25570116989825814, (3, 3) = .0, (3, 4) = .0, (3, 5) = .0, (4, 1) = 3.3148251870684886, (4, 2) = 2.896124015972123, (4, 3) = .9986419139977808, (4, 4) = .0, (4, 5) = .0, (5, 1) = 1.2212245092262748, (5, 2) = 6.019134481287752, (5, 3) = 12.537083329320874, (5, 4) = -.687886036105895, (5, 5) = .0, (6, 1) = 1.2212245092262748, (6, 2) = 6.019134481287752, (6, 3) = 12.537083329320874, (6, 4) = -.687886036105895, (6, 5) = 1.0}, datatype = float[8], order = C_order), Array(1..6, 1..5, {(1, 1) = .0, (1, 2) = .0, (1, 3) = .0, (1, 4) = .0, (1, 5) = .0, (2, 1) = -5.6688, (2, 2) = .0, (2, 3) = .0, (2, 4) = .0, (2, 5) = .0, (3, 1) = -2.4300933568337584, (3, 2) = -.20635991570891224, (3, 3) = .0, (3, 4) = .0, (3, 5) = .0, (4, 1) = -.10735290581452621, (4, 2) = -9.594562251021896, (4, 3) = -20.470286148096154, (4, 4) = .0, (4, 5) = .0, (5, 1) = 7.496443313968615, (5, 2) = -10.246804314641219, (5, 3) = -33.99990352819906, (5, 4) = 11.708908932061595, (5, 5) = .0, (6, 1) = 8.083246795922411, (6, 2) = -7.981132988062785, (6, 3) = -31.52159432874373, (6, 4) = 16.319305431231363, (6, 5) = -6.0588182388340535}, datatype = float[8], order = C_order), Array(1..3, 1..5, {(1, 1) = .0, (1, 2) = .0, (1, 3) = .0, (1, 4) = .0, (1, 5) = .0, (2, 1) = 10.126235083446911, (2, 2) = -7.487995877607633, (2, 3) = -34.800918615557414, (2, 4) = -7.9927717075687275, (2, 5) = 1.0251377232956207, (3, 1) = -.6762803392806898, (3, 2) = 6.087714651678606, (3, 3) = 16.43084320892463, (3, 4) = 24.767225114183653, (3, 5) = -6.5943891257167815}, datatype = float[8], order = C_order)]), ( 9 ) = ([Array(1..3, {(1) = .1, (2) = .1, (3) = .1}, datatype = float[8], order = C_order), Array(1..3, {(1) = .0, (2) = .0, (3) = .0}, datatype = float[8], order = C_order), Array(1..3, {(1) = .0, (2) = .0, (3) = .0}, datatype = float[8], order = C_order), Array(1..3, {(1) = .0, (2) = .0, (3) = .0}, datatype = float[8], order = C_order), Array(1..3, {(1) = .0, (2) = .0, (3) = .0}, datatype = float[8], order = C_order), Array(1..3, 1..3, {(1, 1) = .0, (1, 2) = .0, (1, 3) = .0, (2, 1) = .0, (2, 2) = .0, (2, 3) = .0, (3, 1) = .0, (3, 2) = .0, (3, 3) = .0}, datatype = float[8], order = C_order), Array(1..3, 1..3, {(1, 1) = .0, (1, 2) = .0, (1, 3) = .0, (2, 1) = .0, (2, 2) = .0, (2, 3) = .0, (3, 1) = .0, (3, 2) = .0, (3, 3) = .0}, datatype = float[8], order = C_order), Array(1..3, 1..6, {(1, 1) = .0, (1, 2) = .0, (1, 3) = .0, (1, 4) = .0, (1, 5) = .0, (1, 6) = .0, (2, 1) = .0, (2, 2) = .0, (2, 3) = .0, (2, 4) = .0, (2, 5) = .0, (2, 6) = .0, (3, 1) = .0, (3, 2) = .0, (3, 3) = .0, (3, 4) = .0, (3, 5) = .0, (3, 6) = .0}, datatype = float[8], order = C_order), Array(1..3, {(1) = 0, (2) = 0, (3) = 0}, datatype = integer[8]), Array(1..3, {(1) = .0, (2) = .0, (3) = .0}, datatype = float[8], order = C_order), Array(1..3, {(1) = .0, (2) = .0, (3) = .0}, datatype = float[8], order = C_order), Array(1..3, {(1) = .0, (2) = .0, (3) = .0}, datatype = float[8], order = C_order), Array(1..3, {(1) = .0, (2) = .0, (3) = .0}, datatype = float[8], order = C_order), Array(1..3, {(1) = .0, (2) = .0, (3) = .0}, datatype = float[8], order = C_order)]), ( 8 ) = ([Array(1..3, {(1) = .2342493224, (2) = .1799302829, (3) = .5137385362}, datatype = float[8], order = C_order), Array(1..3, {(1) = .0, (2) = .0, (3) = .0}, datatype = float[8], order = C_order), Array(1..3, {(1) = -.6534884149815456, (2) = .26224111021755525, (3) = 1.1915897756055451}, datatype = float[8], order = C_order), 0, 0]), ( 11 ) = (Array(1..6, 0..3, {(1, 1) = .0, (1, 2) = .0, (1, 3) = .0, (2, 0) = .0, (2, 1) = .0, (2, 2) = .0, (2, 3) = .0, (3, 0) = .0, (3, 1) = .0, (3, 2) = .0, (3, 3) = .0, (4, 0) = .0, (4, 1) = .0, (4, 2) = .0, (4, 3) = .0, (5, 0) = .0, (5, 1) = .0, (5, 2) = .0, (5, 3) = .0, (6, 0) = .0, (6, 1) = .0, (6, 2) = .0, (6, 3) = .0}, datatype = float[8], order = C_order)), ( 10 ) = ([proc (N, X, Y, YP) option `[Y[1] = P(t), Y[2] = Q(t), Y[3] = R(t)]`; YP[1] := -(1/2)*Y[3]-3.46410161513776*X(X)*Y[2]-(1.15470053837925*(-Y(X)^2+2)*X(X)+1.73205080756888*(-2*X(X)^2+Z(X)+1)*Y(X))*Y[1]; YP[2] := (2/3)*Y[3]-1.15470053837925*((1/3)*Y(X)+X(X))*Y[1]+(-1.15470053837925*(-Y(X)^2+2)*X(X)+3.46410161513776*(X(X)^2-(1/2)*Z(X)-2/3)*X(X))*Y[2]; YP[3] := Y[1]*Z(X)+(-1.15470053837925*(-Y(X)^2+2)*X(X)+1.73205080756888*(-2*X(X)^2+Z(X)+4/3)*Y(X))*Y[3]; 0 end proc, -1, 0, 0, 0, 0, 0, 0]), ( 13 ) = (), ( 12 ) = (), ( 15 ) = ("rkf45"), ( 14 ) = ([0, 0]), ( 18 ) = ([]), ( 19 ) = (0), ( 16 ) = ([0, 0, 0, []]), ( 17 ) = ([proc (N, X, Y, YP) option `[Y[1] = P(t), Y[2] = Q(t), Y[3] = R(t)]`; YP[1] := -(1/2)*Y[3]-3.46410161513776*X(X)*Y[2]-(1.15470053837925*(-Y(X)^2+2)*X(X)+1.73205080756888*(-2*X(X)^2+Z(X)+1)*Y(X))*Y[1]; YP[2] := (2/3)*Y[3]-1.15470053837925*((1/3)*Y(X)+X(X))*Y[1]+(-1.15470053837925*(-Y(X)^2+2)*X(X)+3.46410161513776*(X(X)^2-(1/2)*Z(X)-2/3)*X(X))*Y[2]; YP[3] := Y[1]*Z(X)+(-1.15470053837925*(-Y(X)^2+2)*X(X)+1.73205080756888*(-2*X(X)^2+Z(X)+4/3)*Y(X))*Y[3]; 0 end proc, -1, 0, 0, 0, 0, 0, 0]), ( 22 ) = (0), ( 23 ) = (0), ( 20 ) = ([]), ( 21 ) = (0), ( 24 ) = (0) ] )) ] ); _y0 := Array(0..3, {(1) = 0., (2) = .2342493224, (3) = .1799302829}); _vmap := array( 1 .. 3, [( 1 ) = (1), ( 2 ) = (2), ( 3 ) = (3) ] ); _x0 := _dtbl[1][5][5]; _n := _dtbl[1][4][1]; _ne := _dtbl[1][4][3]; _nd := _dtbl[1][4][4]; _nv := _dtbl[1][4][16]; if not type(_xout, 'numeric') then if member(_xout, ["start", "left", "right"]) then if _Env_smart_dsolve_numeric = true or _dtbl[1][4][10] = 1 then if _xout = "left" then if type(_dtbl[2], 'table') then return _dtbl[2][5][1] end if elif _xout = "right" then if type(_dtbl[3], 'table') then return _dtbl[3][5][1] end if end if end if; return _dtbl[1][5][5] elif _xout = "method" then return _dtbl[1][15] elif _xout = "storage" then return evalb(_dtbl[1][4][10] = 1) elif _xout = "leftdata" then if not type(_dtbl[2], 'array') then return NULL else return eval(_dtbl[2]) end if elif _xout = "rightdata" then if not type(_dtbl[3], 'array') then return NULL else return eval(_dtbl[3]) end if elif _xout = "enginedata" then return eval(_dtbl[1]) elif _xout = "enginereset" then _dtbl[2] := evaln(_dtbl[2]); _dtbl[3] := evaln(_dtbl[3]); return NULL elif _xout = "initial" then return procname(_y0[0]) elif _xout = "laxtol" then return _dtbl[`if`(member(_dtbl[4], {2, 3}), _dtbl[4], 1)][5][18] elif _xout = "numfun" then return `if`(member(_dtbl[4], {2, 3}), _dtbl[_dtbl[4]][4][18], 0) elif _xout = "parameters" then return [seq(_y0[_n+_i], _i = 1 .. nops(_pars))] elif _xout = "initial_and_parameters" then return procname(_y0[0]), [seq(_y0[_n+_i], _i = 1 .. nops(_pars))] elif _xout = "last" then if _dtbl[4] <> 2 and _dtbl[4] <> 3 or _x0-_dtbl[_dtbl[4]][5][1] = 0. then error "no information is available on last computed point" else _xout := _dtbl[_dtbl[4]][5][1] end if elif _xout = "function" then if _dtbl[1][4][33]-2. = 0 then return eval(_dtbl[1][10], 1) else return eval(_dtbl[1][10][1], 1) end if elif _xout = "map" then return copy(_vmap) elif type(_xin, `=`) and type(rhs(_xin), 'list') and member(lhs(_xin), {"initial", "parameters", "initial_and_parameters"}) then _ini, _par := [], []; if lhs(_xin) = "initial" then _ini := rhs(_xin) elif lhs(_xin) = "parameters" then _par := rhs(_xin) elif select(type, rhs(_xin), `=`) <> [] then _par, _ini := selectremove(type, rhs(_xin), `=`) elif nops(rhs(_xin)) < nops(_pars)+1 then error "insufficient data for specification of initial and parameters" else _par := rhs(_xin)[-nops(_pars) .. -1]; _ini := rhs(_xin)[1 .. -nops(_pars)-1] end if; _xout := lhs(_xout); if _par <> [] then `dsolve/numeric/process_parameters`(_n, _pars, _par, _y0) end if; if _ini <> [] then `dsolve/numeric/process_initial`(_n-_ne, _ini, _y0, _pars, _vmap) end if; `dsolve/numeric/SC/reinitialize`(_dtbl, _y0, _n, procname, _pars); if _Env_smart_dsolve_numeric = true and type(_y0[0], 'numeric') and _dtbl[1][4][10] <> 1 then procname("right") := _y0[0]; procname("left") := _y0[0] end if; if _xout = "initial" then return [_y0[0], seq(_y0[_vmap[_i]], _i = 1 .. _n-_ne)] elif _xout = "parameters" then return [seq(_y0[_n+_i], _i = 1 .. nops(_pars))] else return [_y0[0], seq(_y0[_vmap[_i]], _i = 1 .. _n-_ne)], [seq(_y0[_n+_i], _i = 1 .. nops(_pars))] end if elif _xin = "eventstop" then if _nv = 0 then error "this solution has no events" end if; _i := _dtbl[4]; if _i <> 2 and _i <> 3 then return 0 end if; if _dtbl[_i][4][10] = 1 and assigned(_dtbl[5-_i]) and _dtbl[_i][4][9] < 100 and 100 <= _dtbl[5-_i][4][9] then _i := 5-_i; _dtbl[4] := _i; _j := round(_dtbl[_i][4][17]); return round(_dtbl[_i][3][1][_j, 1]) elif 100 <= _dtbl[_i][4][9] then _j := round(_dtbl[_i][4][17]); return round(_dtbl[_i][3][1][_j, 1]) else return 0 end if elif _xin = "eventstatus" then if _nv = 0 then error "this solution has no events" end if; _i := [selectremove(proc (a) options operator, arrow; _dtbl[1][3][1][a, 7] = 1 end proc, {seq(_j, _j = 1 .. round(_dtbl[1][3][1][_nv+1, 1]))})]; return ':-enabled' = _i[1], ':-disabled' = _i[2] elif _xin = "eventclear" then if _nv = 0 then error "this solution has no events" end if; _i := _dtbl[4]; if _i <> 2 and _i <> 3 then error "no events to clear" end if; if _dtbl[_i][4][10] = 1 and assigned(_dtbl[5-_i]) and _dtbl[_i][4][9] < 100 and 100 < _dtbl[5-_i][4][9] then _dtbl[4] := 5-_i; _i := 5-_i end if; if _dtbl[_i][4][9] < 100 then error "no events to clear" elif _nv < _dtbl[_i][4][9]-100 then error "event error condition cannot be cleared" else _j := _dtbl[_i][4][9]-100; if irem(round(_dtbl[_i][3][1][_j, 4]), 2) = 1 then error "retriggerable events cannot be cleared" end if; _j := round(_dtbl[_i][3][1][_j, 1]); for _k to _nv do if _dtbl[_i][3][1][_k, 1] = _j then if _dtbl[_i][3][1][_k, 2] = 3 then error "range events cannot be cleared" end if; _dtbl[_i][3][1][_k, 8] := _dtbl[_i][3][1][_nv+1, 8] end if end do; _dtbl[_i][4][17] := 0; _dtbl[_i][4][9] := 0; if _dtbl[1][4][10] = 1 then if _i = 2 then try procname(procname("left")) catch: end try else try procname(procname("right")) catch: end try end if end if end if; return elif type(_xin, `=`) and member(lhs(_xin), {"eventdisable", "eventenable"}) then if _nv = 0 then error "this solution has no events" end if; if type(rhs(_xin), {('list')('posint'), ('set')('posint')}) then _i := {op(rhs(_xin))} elif type(rhs(_xin), 'posint') then _i := {rhs(_xin)} else error "event identifiers must be integers in the range 1..%1", round(_dtbl[1][3][1][_nv+1, 1]) end if; if select(proc (a) options operator, arrow; _nv < a end proc, _i) <> {} then error "event identifiers must be integers in the range 1..%1", round(_dtbl[1][3][1][_nv+1, 1]) end if; _k := {}; for _j to _nv do if member(round(_dtbl[1][3][1][_j, 1]), _i) then _k := `union`(_k, {_j}) end if end do; _i := _k; if lhs(_xin) = "eventdisable" then _dtbl[4] := 0; _j := [evalb(assigned(_dtbl[2]) and member(_dtbl[2][4][17], _i)), evalb(assigned(_dtbl[3]) and member(_dtbl[3][4][17], _i))]; for _k in _i do _dtbl[1][3][1][_k, 7] := 0; if assigned(_dtbl[2]) then _dtbl[2][3][1][_k, 7] := 0 end if; if assigned(_dtbl[3]) then _dtbl[3][3][1][_k, 7] := 0 end if end do; if _j[1] then for _k to _nv+1 do if _k <= _nv and not type(_dtbl[2][3][4][_k, 1], 'undefined') then userinfo(3, {'events', 'eventreset'}, `reinit #2, event code `, _k, ` to defined init `, _dtbl[2][3][4][_k, 1]); _dtbl[2][3][1][_k, 8] := _dtbl[2][3][4][_k, 1] elif _dtbl[2][3][1][_k, 2] = 0 and irem(iquo(round(_dtbl[2][3][1][_k, 4]), 32), 2) = 1 then userinfo(3, {'events', 'eventreset'}, `reinit #2, event code `, _k, ` to rate hysteresis init `, _dtbl[2][5][24]); _dtbl[2][3][1][_k, 8] := _dtbl[2][5][24] elif _dtbl[2][3][1][_k, 2] = 0 and irem(iquo(round(_dtbl[2][3][1][_k, 4]), 2), 2) = 0 then userinfo(3, {'events', 'eventreset'}, `reinit #2, event code `, _k, ` to initial init `, _x0); _dtbl[2][3][1][_k, 8] := _x0 else userinfo(3, {'events', 'eventreset'}, `reinit #2, event code `, _k, ` to fireinitial init `, _x0-1); _dtbl[2][3][1][_k, 8] := _x0-1 end if end do; _dtbl[2][4][17] := 0; _dtbl[2][4][9] := 0; if _dtbl[1][4][10] = 1 then procname(procname("left")) end if end if; if _j[2] then for _k to _nv+1 do if _k <= _nv and not type(_dtbl[3][3][4][_k, 2], 'undefined') then userinfo(3, {'events', 'eventreset'}, `reinit #3, event code `, _k, ` to defined init `, _dtbl[3][3][4][_k, 2]); _dtbl[3][3][1][_k, 8] := _dtbl[3][3][4][_k, 2] elif _dtbl[3][3][1][_k, 2] = 0 and irem(iquo(round(_dtbl[3][3][1][_k, 4]), 32), 2) = 1 then userinfo(3, {'events', 'eventreset'}, `reinit #3, event code `, _k, ` to rate hysteresis init `, _dtbl[3][5][24]); _dtbl[3][3][1][_k, 8] := _dtbl[3][5][24] elif _dtbl[3][3][1][_k, 2] = 0 and irem(iquo(round(_dtbl[3][3][1][_k, 4]), 2), 2) = 0 then userinfo(3, {'events', 'eventreset'}, `reinit #3, event code `, _k, ` to initial init `, _x0); _dtbl[3][3][1][_k, 8] := _x0 else userinfo(3, {'events', 'eventreset'}, `reinit #3, event code `, _k, ` to fireinitial init `, _x0+1); _dtbl[3][3][1][_k, 8] := _x0+1 end if end do; _dtbl[3][4][17] := 0; _dtbl[3][4][9] := 0; if _dtbl[1][4][10] = 1 then procname(procname("right")) end if end if else for _k in _i do _dtbl[1][3][1][_k, 7] := 1 end do; _dtbl[2] := evaln(_dtbl[2]); _dtbl[3] := evaln(_dtbl[3]); _dtbl[4] := 0; if _dtbl[1][4][10] = 1 then if _x0 <= procname("right") then try procname(procname("right")) catch: end try end if; if procname("left") <= _x0 then try procname(procname("left")) catch: end try end if end if end if; return elif type(_xin, `=`) and lhs(_xin) = "eventfired" then if not type(rhs(_xin), 'list') then error "'eventfired' must be specified as a list" end if; if _nv = 0 then error "this solution has no events" end if; if _dtbl[4] <> 2 and _dtbl[4] <> 3 then error "'direction' must be set prior to calling/setting 'eventfired'" end if; _i := _dtbl[4]; _val := NULL; if not assigned(_EnvEventRetriggerWarned) then _EnvEventRetriggerWarned := false end if; for _k in rhs(_xin) do if type(_k, 'integer') then _src := _k elif type(_k, 'integer' = 'anything') and type(evalf(rhs(_k)), 'numeric') then _k := lhs(_k) = evalf[max(Digits, 18)](rhs(_k)); _src := lhs(_k) else error "'eventfired' entry is not valid: %1", _k end if; if _src < 1 or round(_dtbl[1][3][1][_nv+1, 1]) < _src then error "event identifiers must be integers in the range 1..%1", round(_dtbl[1][3][1][_nv+1, 1]) end if; _src := {seq(`if`(_dtbl[1][3][1][_j, 1]-_src = 0., _j, NULL), _j = 1 .. _nv)}; if nops(_src) <> 1 then error "'eventfired' can only be set/queried for root-finding events and time/interval events" end if; _src := _src[1]; if _dtbl[1][3][1][_src, 2] <> 0. and _dtbl[1][3][1][_src, 2]-2. <> 0. then error "'eventfired' can only be set/queried for root-finding events and time/interval events" elif irem(round(_dtbl[1][3][1][_src, 4]), 2) = 1 then if _EnvEventRetriggerWarned = false then WARNING(`'eventfired' has no effect on events that retrigger`) end if; _EnvEventRetriggerWarned := true end if; if _dtbl[_i][3][1][_src, 2] = 0 and irem(iquo(round(_dtbl[_i][3][1][_src, 4]), 32), 2) = 1 then _val := _val, undefined elif type(_dtbl[_i][3][4][_src, _i-1], 'undefined') or _i = 2 and _dtbl[2][3][1][_src, 8] < _dtbl[2][3][4][_src, 1] or _i = 3 and _dtbl[3][3][4][_src, 2] < _dtbl[3][3][1][_src, 8] then _val := _val, _dtbl[_i][3][1][_src, 8] else _val := _val, _dtbl[_i][3][4][_src, _i-1] end if; if type(_k, `=`) then if _dtbl[_i][3][1][_src, 2] = 0 and irem(iquo(round(_dtbl[_i][3][1][_src, 4]), 32), 2) = 1 then error "cannot set event code for a rate hysteresis event" end if; userinfo(3, {'events', 'eventreset'}, `manual set event code `, _src, ` to value `, rhs(_k)); _dtbl[_i][3][1][_src, 8] := rhs(_k); _dtbl[_i][3][4][_src, _i-1] := rhs(_k) end if end do; return [_val] elif type(_xin, `=`) and lhs(_xin) = "direction" then if not member(rhs(_xin), {-1, 1, ':-left', ':-right'}) then error "'direction' must be specified as either '1' or 'right' (positive) or '-1' or 'left' (negative)" end if; _src := `if`(_dtbl[4] = 2, -1, `if`(_dtbl[4] = 3, 1, undefined)); _i := `if`(member(rhs(_xin), {1, ':-right'}), 3, 2); _dtbl[4] := _i; _dtbl[_i] := `dsolve/numeric/SC/IVPdcopy`(_dtbl[1], `if`(assigned(_dtbl[_i]), _dtbl[_i], NULL)); if 0 < _nv then for _j to _nv+1 do if _j <= _nv and not type(_dtbl[_i][3][4][_j, _i-1], 'undefined') then userinfo(3, {'events', 'eventreset'}, `reinit #4, event code `, _j, ` to defined init `, _dtbl[_i][3][4][_j, _i-1]); _dtbl[_i][3][1][_j, 8] := _dtbl[_i][3][4][_j, _i-1] elif _dtbl[_i][3][1][_j, 2] = 0 and irem(iquo(round(_dtbl[_i][3][1][_j, 4]), 32), 2) = 1 then userinfo(3, {'events', 'eventreset'}, `reinit #4, event code `, _j, ` to rate hysteresis init `, _dtbl[_i][5][24]); _dtbl[_i][3][1][_j, 8] := _dtbl[_i][5][24] elif _dtbl[_i][3][1][_j, 2] = 0 and irem(iquo(round(_dtbl[_i][3][1][_j, 4]), 2), 2) = 0 then userinfo(3, {'events', 'eventreset'}, `reinit #4, event code `, _j, ` to initial init `, _x0); _dtbl[_i][3][1][_j, 8] := _x0 else userinfo(3, {'events', 'eventreset'}, `reinit #4, event code `, _j, ` to fireinitial init `, _x0-2*_i+5.0); _dtbl[_i][3][1][_j, 8] := _x0-2*_i+5.0 end if end do end if; return _src elif _xin = "eventcount" then if _dtbl[1][3][1] = 0 or _dtbl[4] <> 2 and _dtbl[4] <> 3 then return 0 else return round(_dtbl[_dtbl[4]][3][1][_nv+1, 12]) end if else return "procname" end if end if; if _xout = _x0 then return [_x0, seq(evalf(_dtbl[1][6][_vmap[_i]]), _i = 1 .. _n-_ne)] end if; _i := `if`(_x0 <= _xout, 3, 2); if _xin = "last" and 0 < _dtbl[_i][4][9] and _dtbl[_i][4][9] < 100 then _dat := eval(_dtbl[_i], 2); _j := _dat[4][20]; return [_dat[11][_j, 0], seq(_dat[11][_j, _vmap[_i]], _i = 1 .. _n-_ne-_nd), seq(_dat[8][1][_vmap[_i]], _i = _n-_ne-_nd+1 .. _n-_ne)] end if; if not type(_dtbl[_i], 'array') then _dtbl[_i] := `dsolve/numeric/SC/IVPdcopy`(_dtbl[1], `if`(assigned(_dtbl[_i]), _dtbl[_i], NULL)); if 0 < _nv then for _j to _nv+1 do if _j <= _nv and not type(_dtbl[_i][3][4][_j, _i-1], 'undefined') then userinfo(3, {'events', 'eventreset'}, `reinit #5, event code `, _j, ` to defined init `, _dtbl[_i][3][4][_j, _i-1]); _dtbl[_i][3][1][_j, 8] := _dtbl[_i][3][4][_j, _i-1] elif _dtbl[_i][3][1][_j, 2] = 0 and irem(iquo(round(_dtbl[_i][3][1][_j, 4]), 32), 2) = 1 then userinfo(3, {'events', 'eventreset'}, `reinit #5, event code `, _j, ` to rate hysteresis init `, _dtbl[_i][5][24]); _dtbl[_i][3][1][_j, 8] := _dtbl[_i][5][24] elif _dtbl[_i][3][1][_j, 2] = 0 and irem(iquo(round(_dtbl[_i][3][1][_j, 4]), 2), 2) = 0 then userinfo(3, {'events', 'eventreset'}, `reinit #5, event code `, _j, ` to initial init `, _x0); _dtbl[_i][3][1][_j, 8] := _x0 else userinfo(3, {'events', 'eventreset'}, `reinit #5, event code `, _j, ` to fireinitial init `, _x0-2*_i+5.0); _dtbl[_i][3][1][_j, 8] := _x0-2*_i+5.0 end if end do end if end if; if _xin <> "last" then if 0 < 0 then if `dsolve/numeric/checkglobals`(op(_dtbl[1][14]), _pars, _n, _y0) then `dsolve/numeric/SC/reinitialize`(_dtbl, _y0, _n, procname, _pars, _i) end if end if; if _dtbl[1][4][7] = 0 then error "parameters must be initialized before solution can be computed" end if end if; _dat := eval(_dtbl[_i], 2); _dtbl[4] := _i; try _src := `dsolve/numeric/SC/IVPrun`(_dat, _xout) catch: userinfo(2, `dsolve/debug`, print(`Exception in solnproc:`, [lastexception][2 .. -1])); error end try; if _src = 0 and 100 < _dat[4][9] then _val := _dat[3][1][_nv+1, 8] else _val := _dat[11][_dat[4][20], 0] end if; if _src <> 0 or _dat[4][9] <= 0 then _dtbl[1][5][1] := _xout else _dtbl[1][5][1] := _val end if; if _i = 3 and _val < _xout then Rounding := -infinity; if _dat[4][9] = 1 then error "cannot evaluate the solution further right of %1, probably a singularity", evalf[8](_val) elif _dat[4][9] = 2 then error "cannot evaluate the solution further right of %1, maxfun limit exceeded (see ?dsolve,maxfun for details)", evalf[8](_val) elif _dat[4][9] = 3 then if _dat[4][25] = 3 then error "cannot evaluate the solution past the initial point, problem may be initially singular or improperly set up" else error "cannot evaluate the solution past the initial point, problem may be complex, initially singular or improperly set up" end if elif _dat[4][9] = 4 then error "cannot evaluate the solution further right of %1, accuracy goal cannot be achieved with specified 'minstep'", evalf[8](_val) elif _dat[4][9] = 5 then error "cannot evaluate the solution further right of %1, too many step failures, tolerances may be too loose for problem", evalf[8](_val) elif _dat[4][9] = 6 then error "cannot evaluate the solution further right of %1, cannot downgrade delay storage for problems with delay derivative order > 1, try increasing delaypts", evalf[8](_val) elif _dat[4][9] = 10 then error "cannot evaluate the solution further right of %1, interrupt requested", evalf[8](_val) elif 100 < _dat[4][9] then if _dat[4][9]-100 = _nv+1 then error "constraint projection failure on event at t=%1", evalf[8](_val) elif _dat[4][9]-100 = _nv+2 then error "index-1 and derivative evaluation failure on event at t=%1", evalf[8](_val) elif _dat[4][9]-100 = _nv+3 then error "maximum number of event iterations reached (%1) at t=%2", round(_dat[3][1][_nv+1, 3]), evalf[8](_val) else if _Env_dsolve_nowarnstop <> true then `dsolve/numeric/warning`(StringTools:-FormatMessage("cannot evaluate the solution further right of %1, event #%2 triggered a halt", evalf[8](_val), round(_dat[3][1][_dat[4][9]-100, 1]))) end if; Rounding := 'nearest'; _xout := _val end if else error "cannot evaluate the solution further right of %1", evalf[8](_val) end if elif _i = 2 and _xout < _val then Rounding := infinity; if _dat[4][9] = 1 then error "cannot evaluate the solution further left of %1, probably a singularity", evalf[8](_val) elif _dat[4][9] = 2 then error "cannot evaluate the solution further left of %1, maxfun limit exceeded (see ?dsolve,maxfun for details)", evalf[8](_val) elif _dat[4][9] = 3 then if _dat[4][25] = 3 then error "cannot evaluate the solution past the initial point, problem may be initially singular or improperly set up" else error "cannot evaluate the solution past the initial point, problem may be complex, initially singular or improperly set up" end if elif _dat[4][9] = 4 then error "cannot evaluate the solution further left of %1, accuracy goal cannot be achieved with specified 'minstep'", evalf[8](_val) elif _dat[4][9] = 5 then error "cannot evaluate the solution further left of %1, too many step failures, tolerances may be too loose for problem", evalf[8](_val) elif _dat[4][9] = 6 then error "cannot evaluate the solution further left of %1, cannot downgrade delay storage for problems with delay derivative order > 1, try increasing delaypts", evalf[8](_val) elif _dat[4][9] = 10 then error "cannot evaluate the solution further right of %1, interrupt requested", evalf[8](_val) elif 100 < _dat[4][9] then if _dat[4][9]-100 = _nv+1 then error "constraint projection failure on event at t=%1", evalf[8](_val) elif _dat[4][9]-100 = _nv+2 then error "index-1 and derivative evaluation failure on event at t=%1", evalf[8](_val) elif _dat[4][9]-100 = _nv+3 then error "maximum number of event iterations reached (%1) at t=%2", round(_dat[3][1][_nv+1, 3]), evalf[8](_val) else if _Env_dsolve_nowarnstop <> true then `dsolve/numeric/warning`(StringTools:-FormatMessage("cannot evaluate the solution further left of %1, event #%2 triggered a halt", evalf[8](_val), round(_dat[3][1][_dat[4][9]-100, 1]))) end if; Rounding := 'nearest'; _xout := _val end if else error "cannot evaluate the solution further left of %1", evalf[8](_val) end if end if; if _EnvInFsolve = true then _dig := _dat[4][26]; _dat[4][26] := _EnvDSNumericSaveDigits; _Env_dsolve_SC_native := true; if _dat[4][25] = 1 then _i := 1; _dat[4][25] := 2 else _i := _dat[4][25] end if; _val := `dsolve/numeric/SC/IVPval`(_dat, _xout, _src); _dat[4][25] := _i; _dat[4][26] := _dig; [_xout, seq(_val[_vmap[_i]], _i = 1 .. _n-_ne)] else Digits := _dat[4][26]; _val := `dsolve/numeric/SC/IVPval`(eval(_dat, 2), _xout, _src); [_xout, seq(_val[_vmap[_i]], _i = 1 .. _n-_ne)] end if end proc, (2) = Array(0..0, {}), (3) = [t, P(t), Q(t), R(t)], (4) = []}); _vars := _dat[3]; _pars := map(rhs, _dat[4]); _n := nops(_vars)-1; _solnproc := _dat[1]; if not type(_xout, 'numeric') then if member(x_rkf45, ["start", 'start', "method", 'method', "left", 'left', "right", 'right', "leftdata", "rightdata", "enginedata", "eventstop", 'eventstop', "eventclear", 'eventclear', "eventstatus", 'eventstatus', "eventcount", 'eventcount', "laxtol", 'laxtol', "numfun", 'numfun', NULL]) then _res := _solnproc(convert(x_rkf45, 'string')); if 1 < nops([_res]) then return _res elif type(_res, 'array') then return eval(_res, 1) elif _res <> "procname" then return _res end if elif member(x_rkf45, ["last", 'last', "initial", 'initial', "parameters", 'parameters', "initial_and_parameters", 'initial_and_parameters', NULL]) then _xout := convert(x_rkf45, 'string'); _res := _solnproc(_xout); if _xout = "parameters" then return [seq(_pars[_i] = _res[_i], _i = 1 .. nops(_pars))] elif _xout = "initial_and_parameters" then return [seq(_vars[_i+1] = [_res][1][_i+1], _i = 0 .. _n), seq(_pars[_i] = [_res][2][_i], _i = 1 .. nops(_pars))] else return [seq(_vars[_i+1] = _res[_i+1], _i = 0 .. _n)] end if elif type(_xout, `=`) and member(lhs(_xout), ["initial", 'initial', "parameters", 'parameters', "initial_and_parameters", 'initial_and_parameters', NULL]) then _xout := convert(lhs(x_rkf45), 'string') = rhs(x_rkf45); if type(rhs(_xout), 'list') then _res := _solnproc(_xout) else error "initial and/or parameter values must be specified in a list" end if; if lhs(_xout) = "initial" then return [seq(_vars[_i+1] = _res[_i+1], _i = 0 .. _n)] elif lhs(_xout) = "parameters" then return [seq(_pars[_i] = _res[_i], _i = 1 .. nops(_pars))] else return [seq(_vars[_i+1] = [_res][1][_i+1], _i = 0 .. _n), seq(_pars[_i] = [_res][2][_i], _i = 1 .. nops(_pars))] end if elif type(_xout, `=`) and member(lhs(_xout), ["eventdisable", 'eventdisable', "eventenable", 'eventenable', "eventfired", 'eventfired', "direction", 'direction', NULL]) then return _solnproc(convert(lhs(x_rkf45), 'string') = rhs(x_rkf45)) elif _xout = "solnprocedure" then return eval(_solnproc) elif _xout = "sysvars" then return _vars end if; if procname <> unknown then return ('procname')(x_rkf45) else _ndsol; _ndsol := pointto(_dat[2][0]); return ('_ndsol')(x_rkf45) end if end if; try _res := _solnproc(_xout); [seq(_vars[_i+1] = _res[_i+1], _i = 0 .. _n)] catch: error end try end proc

 (5)
> 

# All in a row

soln := dsolve({eq1, eq2, eq3, eq4, eq5, eq6, ics, ICS}, {P(t), Q(t), R(t), x(t), y(t), z(t)}, numeric, start = 0, range = 0 .. 1)

proc (x_rkf45) local _res, _dat, _vars, _solnproc, _xout, _ndsol, _pars, _n, _i; option `Copyright (c) 2000 by Waterloo Maple Inc. All rights reserved.`; if 1 < nargs then error "invalid input: too many arguments" end if; _EnvDSNumericSaveDigits := Digits; Digits := 15; if _EnvInFsolve = true then _xout := evalf[_EnvDSNumericSaveDigits](x_rkf45) else _xout := evalf(x_rkf45) end if; _dat := Array(1..4, {(1) = proc (_xin) local _xout, _dtbl, _dat, _vmap, _x0, _y0, _val, _dig, _n, _ne, _nd, _nv, _pars, _ini, _par, _i, _j, _k, _src; option `Copyright (c) 2002 by Waterloo Maple Inc. All rights reserved.`; table( [( "left" ) = 0., ( "complex" ) = false, ( "right" ) = 1. ] ) _xout := _xin; _pars := []; _dtbl := array( 1 .. 4, [( 1 ) = (array( 1 .. 24, [( 1 ) = (datatype = float[8], order = C_order, storage = rectangular), ( 2 ) = (datatype = float[8], order = C_order, storage = rectangular), ( 3 ) = ([0, 0, 0, Array(1..0, 1..2, {}, datatype = float[8], order = C_order)]), ( 4 ) = (Array(1..54, {(1) = 6, (2) = 6, (3) = 0, (4) = 0, (5) = 0, (6) = 0, (7) = 1, (8) = 0, (9) = 0, (10) = 1, (11) = 0, (12) = 0, (13) = 0, (14) = 0, (15) = 0, (16) = 0, (17) = 0, (18) = 1, (19) = 30000, (20) = 0, (21) = 0, (22) = 1, (23) = 4, (24) = 0, (25) = 1, (26) = 15, (27) = 1, (28) = 0, (29) = 1, (30) = 3, (31) = 3, (32) = 0, (33) = 1, (34) = 0, (35) = 0, (36) = 0, (37) = 0, (38) = 0, (39) = 0, (40) = 0, (41) = 0, (42) = 0, (43) = 1, (44) = 0, (45) = 0, (46) = 0, (47) = 0, (48) = 0, (49) = 0, (50) = 50, (51) = 1, (52) = 0, (53) = 0, (54) = 0}, datatype = integer[8])), ( 5 ) = (Array(1..28, {(1) = 1.0, (2) = 0.10e-5, (3) = .0, (4) = 0.500001e-14, (5) = .0, (6) = 0.2581800197474211e-1, (7) = .0, (8) = 0.10e-5, (9) = .0, (10) = .0, (11) = .0, (12) = .0, (13) = 1.0, (14) = .0, (15) = .49999999999999, (16) = .0, (17) = 1.0, (18) = 1.0, (19) = .0, (20) = .0, (21) = 1.0, (22) = 1.0, (23) = .0, (24) = .0, (25) = 0.10e-14, (26) = .0, (27) = .0, (28) = .0}, datatype = float[8], order = C_order)), ( 6 ) = (Array(1..6, {(1) = .2342493224, (2) = .1799302829, (3) = .5137385362, (4) = -0.1e-1, (5) = .99, (6) = 0.1e-1}, datatype = float[8], order = C_order)), ( 7 ) = ([Array(1..4, 1..7, {(1, 1) = .0, (1, 2) = .203125, (1, 3) = .3046875, (1, 4) = .75, (1, 5) = .8125, (1, 6) = .40625, (1, 7) = .8125, (2, 1) = 0.6378173828125e-1, (2, 2) = .0, (2, 3) = .279296875, (2, 4) = .27237892150878906, (2, 5) = -0.9686851501464844e-1, (2, 6) = 0.1956939697265625e-1, (2, 7) = .5381584167480469, (3, 1) = 0.31890869140625e-1, (3, 2) = .0, (3, 3) = -.34375, (3, 4) = -.335235595703125, (3, 5) = .2296142578125, (3, 6) = .41748046875, (3, 7) = 11.480712890625, (4, 1) = 0.9710520505905151e-1, (4, 2) = .0, (4, 3) = .40350341796875, (4, 4) = 0.20297467708587646e-1, (4, 5) = -0.6054282188415527e-2, (4, 6) = -0.4770040512084961e-1, (4, 7) = .77858567237854}, datatype = float[8], order = C_order), Array(1..6, 1..6, {(1, 1) = .0, (1, 2) = .0, (1, 3) = .0, (1, 4) = .0, (1, 5) = .0, (1, 6) = 1.0, (2, 1) = .25, (2, 2) = .0, (2, 3) = .0, (2, 4) = .0, (2, 5) = .0, (2, 6) = 1.0, (3, 1) = .1875, (3, 2) = .5625, (3, 3) = .0, (3, 4) = .0, (3, 5) = .0, (3, 6) = 2.0, (4, 1) = .23583984375, (4, 2) = -.87890625, (4, 3) = .890625, (4, 4) = .0, (4, 5) = .0, (4, 6) = .2681884765625, (5, 1) = .1272735595703125, (5, 2) = -.5009765625, (5, 3) = .44921875, (5, 4) = -0.128936767578125e-1, (5, 5) = .0, (5, 6) = 0.626220703125e-1, (6, 1) = -0.927734375e-1, (6, 2) = .626220703125, (6, 3) = -.4326171875, (6, 4) = .1418304443359375, (6, 5) = -0.861053466796875e-1, (6, 6) = .3131103515625}, datatype = float[8], order = C_order), Array(1..6, {(1) = .0, (2) = .386, (3) = .21, (4) = .63, (5) = 1.0, (6) = 1.0}, datatype = float[8], order = C_order), Array(1..6, {(1) = .25, (2) = -.1043, (3) = .1035, (4) = -0.362e-1, (5) = .0, (6) = .0}, datatype = float[8], order = C_order), Array(1..6, 1..5, {(1, 1) = .0, (1, 2) = .0, (1, 3) = .0, (1, 4) = .0, (1, 5) = .0, (2, 1) = 1.544, (2, 2) = .0, (2, 3) = .0, (2, 4) = .0, (2, 5) = .0, (3, 1) = .9466785280815533, (3, 2) = .25570116989825814, (3, 3) = .0, (3, 4) = .0, (3, 5) = .0, (4, 1) = 3.3148251870684886, (4, 2) = 2.896124015972123, (4, 3) = .9986419139977808, (4, 4) = .0, (4, 5) = .0, (5, 1) = 1.2212245092262748, (5, 2) = 6.019134481287752, (5, 3) = 12.537083329320874, (5, 4) = -.687886036105895, (5, 5) = .0, (6, 1) = 1.2212245092262748, (6, 2) = 6.019134481287752, (6, 3) = 12.537083329320874, (6, 4) = -.687886036105895, (6, 5) = 1.0}, datatype = float[8], order = C_order), Array(1..6, 1..5, {(1, 1) = .0, (1, 2) = .0, (1, 3) = .0, (1, 4) = .0, (1, 5) = .0, (2, 1) = -5.6688, (2, 2) = .0, (2, 3) = .0, (2, 4) = .0, (2, 5) = .0, (3, 1) = -2.4300933568337584, (3, 2) = -.20635991570891224, (3, 3) = .0, (3, 4) = .0, (3, 5) = .0, (4, 1) = -.10735290581452621, (4, 2) = -9.594562251021896, (4, 3) = -20.470286148096154, (4, 4) = .0, (4, 5) = .0, (5, 1) = 7.496443313968615, (5, 2) = -10.246804314641219, (5, 3) = -33.99990352819906, (5, 4) = 11.708908932061595, (5, 5) = .0, (6, 1) = 8.083246795922411, (6, 2) = -7.981132988062785, (6, 3) = -31.52159432874373, (6, 4) = 16.319305431231363, (6, 5) = -6.0588182388340535}, datatype = float[8], order = C_order), Array(1..3, 1..5, {(1, 1) = .0, (1, 2) = .0, (1, 3) = .0, (1, 4) = .0, (1, 5) = .0, (2, 1) = 10.126235083446911, (2, 2) = -7.487995877607633, (2, 3) = -34.800918615557414, (2, 4) = -7.9927717075687275, (2, 5) = 1.0251377232956207, (3, 1) = -.6762803392806898, (3, 2) = 6.087714651678606, (3, 3) = 16.43084320892463, (3, 4) = 24.767225114183653, (3, 5) = -6.5943891257167815}, datatype = float[8], order = C_order)]), ( 9 ) = ([Array(1..6, {(1) = .1, (2) = .1, (3) = .1, (4) = .1, (5) = .1, (6) = .1}, datatype = float[8], order = C_order), Array(1..6, {(1) = -.6469510038574611, (2) = 1.7084449862521758, (3) = 5.619529600793772, (4) = 0.10581330059387134e-2, (5) = .9969522087362933, (6) = 0.16912331444129752e-2}, datatype = float[8], order = C_order), Array(1..6, {(1) = .0, (2) = .0, (3) = .0, (4) = .0, (5) = .0, (6) = .0}, datatype = float[8], order = C_order), Array(1..6, {(1) = -1.8174564552465413, (2) = 4.30366785167187, (3) = 13.960107782236397, (4) = 0.14695938777494755e-2, (5) = 0.33726538182153962e-2, (6) = -0.2756164352142944e-2}, datatype = float[8], order = C_order), Array(1..6, {(1) = .0, (2) = .0, (3) = .0, (4) = .0, (5) = .0, (6) = .0}, datatype = float[8], order = C_order), Array(1..6, 1..6, {(1, 1) = .0, (1, 2) = .0, (1, 3) = .0, (1, 4) = .0, (1, 5) = .0, (1, 6) = .0, (2, 1) = .0, (2, 2) = .0, (2, 3) = .0, (2, 4) = .0, (2, 5) = .0, (2, 6) = .0, (3, 1) = .0, (3, 2) = .0, (3, 3) = .0, (3, 4) = .0, (3, 5) = .0, (3, 6) = .0, (4, 1) = .0, (4, 2) = .0, (4, 3) = .0, (4, 4) = .0, (4, 5) = .0, (4, 6) = .0, (5, 1) = .0, (5, 2) = .0, (5, 3) = .0, (5, 4) = .0, (5, 5) = .0, (5, 6) = .0, (6, 1) = .0, (6, 2) = .0, (6, 3) = .0, (6, 4) = .0, (6, 5) = .0, (6, 6) = .0}, datatype = float[8], order = C_order), Array(1..6, 1..6, {(1, 1) = .0, (1, 2) = .0, (1, 3) = .0, (1, 4) = .0, (1, 5) = .0, (1, 6) = .0, (2, 1) = .0, (2, 2) = .0, (2, 3) = .0, (2, 4) = .0, (2, 5) = .0, (2, 6) = .0, (3, 1) = .0, (3, 2) = .0, (3, 3) = .0, (3, 4) = .0, (3, 5) = .0, (3, 6) = .0, (4, 1) = .0, (4, 2) = .0, (4, 3) = .0, (4, 4) = .0, (4, 5) = .0, (4, 6) = .0, (5, 1) = .0, (5, 2) = .0, (5, 3) = .0, (5, 4) = .0, (5, 5) = .0, (5, 6) = .0, (6, 1) = .0, (6, 2) = .0, (6, 3) = .0, (6, 4) = .0, (6, 5) = .0, (6, 6) = .0}, datatype = float[8], order = C_order), Array(1..6, 1..6, {(1, 1) = -1.8174564552465413, (1, 2) = -1.5961798310015682, (1, 3) = -1.6320984766762368, (1, 4) = -1.794276191365887, (1, 5) = -1.8196239216963577, (1, 6) = -1.666791376880401, (2, 1) = 4.30366785167187, (2, 2) = 3.7285529660791323, (2, 3) = 3.8220571587930636, (2, 4) = 4.2436470774636215, (2, 5) = 4.309254368325473, (2, 6) = 3.9125540306811586, (3, 1) = 13.960107782236397, (3, 2) = 12.107313789961491, (3, 3) = 12.408495351207202, (3, 4) = 13.766676071871467, (3, 5) = 13.978102815543973, (3, 6) = 12.699960113697669, (4, 1) = 0.14695938777494755e-2, (4, 2) = 0.1890417736254646e-2, (4, 3) = 0.18192368255635083e-2, (4, 4) = 0.15081516037365184e-2, (4, 5) = 0.1464218771815906e-2, (4, 6) = 0.17469485986022493e-2, (5, 1) = 0.33726538182153962e-2, (5, 2) = 0.3616834102103271e-2, (5, 3) = 0.3575806202298267e-2, (5, 4) = 0.3396699888512361e-2, (5, 5) = 0.33718646416661044e-2, (5, 6) = 0.3534391454307867e-2, (6, 1) = -0.2756164352142944e-2, (6, 2) = -0.30603548985099053e-2, (6, 3) = -0.3009093417590724e-2, (6, 4) = -0.2785171309676129e-2, (6, 5) = -0.27538569589072135e-2, (6, 6) = -0.2957201051683882e-2}, datatype = float[8], order = C_order), Array(1..6, {(1) = 0, (2) = 0, (3) = 0, (4) = 0, (5) = 0, (6) = 0}, datatype = integer[8]), Array(1..6, {(1) = -.58658912992647, (2) = 1.567090440655536, (3) = 5.160614737970771, (4) = 0.9908058161734378e-3, (5) = .9968199579900292, (6) = 0.18025382964810283e-2}, datatype = float[8], order = C_order), Array(1..6, {(1) = -.7043836227107904, (2) = 1.8439931145266863, (3) = 6.059328015409783, (4) = 0.11097276577890164e-2, (5) = .997064605805955, (6) = 0.15985092092862313e-2}, datatype = float[8], order = C_order), Array(1..6, {(1) = 0.54780767477780046e-7, (2) = 0.1349491964841576e-6, (3) = 0.4366842834002682e-6, (4) = 0.6213567312583768e-10, (5) = 0.6205702618444775e-11, (6) = 0.24046496564789077e-10}, datatype = float[8], order = C_order), Array(1..6, {(1) = .2342493224, (2) = .1799302829, (3) = .5137385362, (4) = -0.1e-1, (5) = .99, (6) = 0.1e-1}, datatype = float[8], order = C_order), Array(1..6, {(1) = -1.7035037265233508, (2) = 4.008091885160755, (3) = 13.007711580584129, (4) = 0.16741346853412775e-2, (5) = 0.34931094516621056e-2, (6) = -0.290522356073048e-2}, datatype = float[8], order = C_order)]), ( 8 ) = ([Array(1..6, {(1) = .2342493224, (2) = .1799302829, (3) = .5137385362, (4) = -0.1e-1, (5) = .99, (6) = 0.1e-1}, datatype = float[8], order = C_order), Array(1..6, {(1) = .0, (2) = .0, (3) = .0, (4) = .0, (5) = .0, (6) = .0}, datatype = float[8], order = C_order), Array(1..6, {(1) = -.6534884149815456, (2) = .26224111021755525, (3) = 1.1915897756055451, (4) = 0.28547973179495756e-1, (5) = 0.1142722829694645e-1, (6) = -0.1697010265031224e-1}, datatype = float[8], order = C_order), 0, 0]), ( 11 ) = (Array(1..6, 0..6, {(1, 1) = .0, (1, 2) = .0, (1, 3) = .0, (1, 4) = .0, (1, 5) = .0, (1, 6) = .0, (2, 0) = .0, (2, 1) = .0, (2, 2) = .0, (2, 3) = .0, (2, 4) = .0, (2, 5) = .0, (2, 6) = .0, (3, 0) = .0, (3, 1) = .0, (3, 2) = .0, (3, 3) = .0, (3, 4) = .0, (3, 5) = .0, (3, 6) = .0, (4, 0) = .0, (4, 1) = .0, (4, 2) = .0, (4, 3) = .0, (4, 4) = .0, (4, 5) = .0, (4, 6) = .0, (5, 0) = .0, (5, 1) = .0, (5, 2) = .0, (5, 3) = .0, (5, 4) = .0, (5, 5) = .0, (5, 6) = .0, (6, 0) = .0, (6, 1) = .0, (6, 2) = .0, (6, 3) = .0, (6, 4) = .0, (6, 5) = .0, (6, 6) = .0}, datatype = float[8], order = C_order)), ( 10 ) = ([proc (N, X, Y, YP) option `[Y[1] = P(t), Y[2] = Q(t), Y[3] = R(t), Y[4] = x(t), Y[5] = y(t), Y[6] = z(t)]`; YP[1] := -(1/2)*Y[3]-3.46410161513776*X(X)*Y[2]-(1.15470053837925*(-Y(X)^2+2)*X(X)+1.73205080756888*(-2*X(X)^2+Z(X)+1)*Y(X))*Y[1]; YP[2] := (2/3)*Y[3]-1.15470053837925*((1/3)*Y(X)+X(X))*Y[1]+(-1.15470053837925*(-Y(X)^2+2)*X(X)+3.46410161513776*(X(X)^2-(1/2)*Z(X)-2/3)*X(X))*Y[2]; YP[3] := Y[1]*Z(X)+(-1.15470053837925*(-Y(X)^2+2)*X(X)+1.73205080756888*(-2*X(X)^2+Z(X)+4/3)*Y(X))*Y[3]; YP[4] := 1.73205080756888*Y[4]^3*Y[5]+.288675134594813*(2*Y[5]^2-2)*Y[4]^2+.866025403784440*Y[5]*(Y[6]-2)*Y[4]-.577350269189627*Y[5]^2+.577350269189627; YP[5] := .288675134594813*(Y[5]-1)*(Y[5]+1)*(6*Y[4]^2+2*Y[4]*Y[5]+3*Y[6]-2); YP[6] := .577350269189627*Y[6]*(6*Y[4]^2*Y[5]+2*Y[4]*Y[5]^2+3*Y[5]*Y[6]-2*Y[4]-3*Y[5]); 0 end proc, -1, 0, 0, 0, 0, 0, 0]), ( 13 ) = (), ( 12 ) = (), ( 15 ) = ("rkf45"), ( 14 ) = ([0, 0]), ( 18 ) = ([]), ( 19 ) = (0), ( 16 ) = ([0, 0, 0, []]), ( 17 ) = ([proc (N, X, Y, YP) option `[Y[1] = P(t), Y[2] = Q(t), Y[3] = R(t), Y[4] = x(t), Y[5] = y(t), Y[6] = z(t)]`; YP[1] := -(1/2)*Y[3]-3.46410161513776*X(X)*Y[2]-(1.15470053837925*(-Y(X)^2+2)*X(X)+1.73205080756888*(-2*X(X)^2+Z(X)+1)*Y(X))*Y[1]; YP[2] := (2/3)*Y[3]-1.15470053837925*((1/3)*Y(X)+X(X))*Y[1]+(-1.15470053837925*(-Y(X)^2+2)*X(X)+3.46410161513776*(X(X)^2-(1/2)*Z(X)-2/3)*X(X))*Y[2]; YP[3] := Y[1]*Z(X)+(-1.15470053837925*(-Y(X)^2+2)*X(X)+1.73205080756888*(-2*X(X)^2+Z(X)+4/3)*Y(X))*Y[3]; YP[4] := 1.73205080756888*Y[4]^3*Y[5]+.288675134594813*(2*Y[5]^2-2)*Y[4]^2+.866025403784440*Y[5]*(Y[6]-2)*Y[4]-.577350269189627*Y[5]^2+.577350269189627; YP[5] := .288675134594813*(Y[5]-1)*(Y[5]+1)*(6*Y[4]^2+2*Y[4]*Y[5]+3*Y[6]-2); YP[6] := .577350269189627*Y[6]*(6*Y[4]^2*Y[5]+2*Y[4]*Y[5]^2+3*Y[5]*Y[6]-2*Y[4]-3*Y[5]); 0 end proc, -1, 0, 0, 0, 0, 0, 0]), ( 22 ) = (0), ( 23 ) = (0), ( 20 ) = ([]), ( 21 ) = (0), ( 24 ) = (0) ] )), ( 3 ) = (array( 1 .. 24, [( 1 ) = (datatype = float[8], order = C_order, storage = rectangular), ( 2 ) = (datatype = float[8], order = C_order, storage = rectangular), ( 3 ) = ([0, 0, 0, Array(1..0, 1..2, {}, datatype = float[8], order = C_order)]), ( 4 ) = (Array(1..54, {(1) = 6, (2) = 6, (3) = 0, (4) = 0, (5) = 0, (6) = 0, (7) = 1, (8) = 1, (9) = 0, (10) = 1, (11) = 57, (12) = 57, (13) = 0, (14) = 0, (15) = 0, (16) = 0, (17) = 0, (18) = 99, (19) = 30000, (20) = 5, (21) = 0, (22) = 1, (23) = 4, (24) = 0, (25) = 1, (26) = 15, (27) = 1, (28) = 0, (29) = 1, (30) = 3, (31) = 3, (32) = 0, (33) = 1, (34) = 0, (35) = 0, (36) = 0, (37) = 0, (38) = 0, (39) = 0, (40) = 0, (41) = 0, (42) = 0, (43) = 1, (44) = 0, (45) = 0, (46) = 0, (47) = 0, (48) = 0, (49) = 0, (50) = 50, (51) = 1, (52) = 0, (53) = 0, (54) = 0}, datatype = integer[8])), ( 5 ) = (Array(1..28, {(1) = 1.0, (2) = 0.10e-5, (3) = 0.8170088713226942e-1, (4) = 0.500001e-14, (5) = .0, (6) = 0.2581800197474211e-1, (7) = .0, (8) = 0.10e-5, (9) = .0, (10) = .0, (11) = .0, (12) = .0, (13) = 1.0, (14) = .0, (15) = .49999999999999, (16) = .0, (17) = 1.0, (18) = 1.0, (19) = .0, (20) = .0, (21) = 1.0, (22) = 1.0, (23) = .0, (24) = .0, (25) = 0.10e-14, (26) = .0, (27) = .0, (28) = .0}, datatype = float[8], order = C_order)), ( 6 ) = (Array(1..6, {(1) = .2342493224, (2) = .1799302829, (3) = .5137385362, (4) = -0.1e-1, (5) = .99, (6) = 0.1e-1}, datatype = float[8], order = C_order)), ( 7 ) = ([Array(1..4, 1..7, {(1, 1) = .0, (1, 2) = .203125, (1, 3) = .3046875, (1, 4) = .75, (1, 5) = .8125, (1, 6) = .40625, (1, 7) = .8125, (2, 1) = 0.6378173828125e-1, (2, 2) = .0, (2, 3) = .279296875, (2, 4) = .27237892150878906, (2, 5) = -0.9686851501464844e-1, (2, 6) = 0.1956939697265625e-1, (2, 7) = .5381584167480469, (3, 1) = 0.31890869140625e-1, (3, 2) = .0, (3, 3) = -.34375, (3, 4) = -.335235595703125, (3, 5) = .2296142578125, (3, 6) = .41748046875, (3, 7) = 11.480712890625, (4, 1) = 0.9710520505905151e-1, (4, 2) = .0, (4, 3) = .40350341796875, (4, 4) = 0.20297467708587646e-1, (4, 5) = -0.6054282188415527e-2, (4, 6) = -0.4770040512084961e-1, (4, 7) = .77858567237854}, datatype = float[8], order = C_order), Array(1..6, 1..6, {(1, 1) = .0, (1, 2) = .0, (1, 3) = .0, (1, 4) = .0, (1, 5) = .0, (1, 6) = 1.0, (2, 1) = .25, (2, 2) = .0, (2, 3) = .0, (2, 4) = .0, (2, 5) = .0, (2, 6) = 1.0, (3, 1) = .1875, (3, 2) = .5625, (3, 3) = .0, (3, 4) = .0, (3, 5) = .0, (3, 6) = 2.0, (4, 1) = .23583984375, (4, 2) = -.87890625, (4, 3) = .890625, (4, 4) = .0, (4, 5) = .0, (4, 6) = .2681884765625, (5, 1) = .1272735595703125, (5, 2) = -.5009765625, (5, 3) = .44921875, (5, 4) = -0.128936767578125e-1, (5, 5) = .0, (5, 6) = 0.626220703125e-1, (6, 1) = -0.927734375e-1, (6, 2) = .626220703125, (6, 3) = -.4326171875, (6, 4) = .1418304443359375, (6, 5) = -0.861053466796875e-1, (6, 6) = .3131103515625}, datatype = float[8], order = C_order), Array(1..6, {(1) = .0, (2) = .386, (3) = .21, (4) = .63, (5) = 1.0, (6) = 1.0}, datatype = float[8], order = C_order), Array(1..6, {(1) = .25, (2) = -.1043, (3) = .1035, (4) = -0.362e-1, (5) = .0, (6) = .0}, datatype = float[8], order = C_order), Array(1..6, 1..5, {(1, 1) = .0, (1, 2) = .0, (1, 3) = .0, (1, 4) = .0, (1, 5) = .0, (2, 1) = 1.544, (2, 2) = .0, (2, 3) = .0, (2, 4) = .0, (2, 5) = .0, (3, 1) = .9466785280815533, (3, 2) = .25570116989825814, (3, 3) = .0, (3, 4) = .0, (3, 5) = .0, (4, 1) = 3.3148251870684886, (4, 2) = 2.896124015972123, (4, 3) = .9986419139977808, (4, 4) = .0, (4, 5) = .0, (5, 1) = 1.2212245092262748, (5, 2) = 6.019134481287752, (5, 3) = 12.537083329320874, (5, 4) = -.687886036105895, (5, 5) = .0, (6, 1) = 1.2212245092262748, (6, 2) = 6.019134481287752, (6, 3) = 12.537083329320874, (6, 4) = -.687886036105895, (6, 5) = 1.0}, datatype = float[8], order = C_order), Array(1..6, 1..5, {(1, 1) = .0, (1, 2) = .0, (1, 3) = .0, (1, 4) = .0, (1, 5) = .0, (2, 1) = -5.6688, (2, 2) = .0, (2, 3) = .0, (2, 4) = .0, (2, 5) = .0, (3, 1) = -2.4300933568337584, (3, 2) = -.20635991570891224, (3, 3) = .0, (3, 4) = .0, (3, 5) = .0, (4, 1) = -.10735290581452621, (4, 2) = -9.594562251021896, (4, 3) = -20.470286148096154, (4, 4) = .0, (4, 5) = .0, (5, 1) = 7.496443313968615, (5, 2) = -10.246804314641219, (5, 3) = -33.99990352819906, (5, 4) = 11.708908932061595, (5, 5) = .0, (6, 1) = 8.083246795922411, (6, 2) = -7.981132988062785, (6, 3) = -31.52159432874373, (6, 4) = 16.319305431231363, (6, 5) = -6.0588182388340535}, datatype = float[8], order = C_order), Array(1..3, 1..5, {(1, 1) = .0, (1, 2) = .0, (1, 3) = .0, (1, 4) = .0, (1, 5) = .0, (2, 1) = 10.126235083446911, (2, 2) = -7.487995877607633, (2, 3) = -34.800918615557414, (2, 4) = -7.9927717075687275, (2, 5) = 1.0251377232956207, (3, 1) = -.6762803392806898, (3, 2) = 6.087714651678606, (3, 3) = 16.43084320892463, (3, 4) = 24.767225114183653, (3, 5) = -6.5943891257167815}, datatype = float[8], order = C_order)]), ( 9 ) = ([Array(1..6, {(1) = .1, (2) = .1, (3) = .1, (4) = .1, (5) = .1, (6) = .1}, datatype = float[8], order = C_order), Array(1..6, {(1) = -.6469510038574611, (2) = 1.7084449862521758, (3) = 5.619529600793772, (4) = 0.10581330059387134e-2, (5) = .9969522087362933, (6) = 0.16912331444129752e-2}, datatype = float[8], order = C_order), Array(1..6, {(1) = .0, (2) = .0, (3) = .0, (4) = .0, (5) = .0, (6) = .0}, datatype = float[8], order = C_order), Array(1..6, {(1) = -1.8174564552465413, (2) = 4.30366785167187, (3) = 13.960107782236397, (4) = 0.14695938777494755e-2, (5) = 0.33726538182153962e-2, (6) = -0.2756164352142944e-2}, datatype = float[8], order = C_order), Array(1..6, {(1) = .0, (2) = .0, (3) = .0, (4) = .0, (5) = .0, (6) = .0}, datatype = float[8], order = C_order), Array(1..6, 1..6, {(1, 1) = .0, (1, 2) = .0, (1, 3) = .0, (1, 4) = .0, (1, 5) = .0, (1, 6) = .0, (2, 1) = .0, (2, 2) = .0, (2, 3) = .0, (2, 4) = .0, (2, 5) = .0, (2, 6) = .0, (3, 1) = .0, (3, 2) = .0, (3, 3) = .0, (3, 4) = .0, (3, 5) = .0, (3, 6) = .0, (4, 1) = .0, (4, 2) = .0, (4, 3) = .0, (4, 4) = .0, (4, 5) = .0, (4, 6) = .0, (5, 1) = .0, (5, 2) = .0, (5, 3) = .0, (5, 4) = .0, (5, 5) = .0, (5, 6) = .0, (6, 1) = .0, (6, 2) = .0, (6, 3) = .0, (6, 4) = .0, (6, 5) = .0, (6, 6) = .0}, datatype = float[8], order = C_order), Array(1..6, 1..6, {(1, 1) = .0, (1, 2) = .0, (1, 3) = .0, (1, 4) = .0, (1, 5) = .0, (1, 6) = .0, (2, 1) = .0, (2, 2) = .0, (2, 3) = .0, (2, 4) = .0, (2, 5) = .0, (2, 6) = .0, (3, 1) = .0, (3, 2) = .0, (3, 3) = .0, (3, 4) = .0, (3, 5) = .0, (3, 6) = .0, (4, 1) = .0, (4, 2) = .0, (4, 3) = .0, (4, 4) = .0, (4, 5) = .0, (4, 6) = .0, (5, 1) = .0, (5, 2) = .0, (5, 3) = .0, (5, 4) = .0, (5, 5) = .0, (5, 6) = .0, (6, 1) = .0, (6, 2) = .0, (6, 3) = .0, (6, 4) = .0, (6, 5) = .0, (6, 6) = .0}, datatype = float[8], order = C_order), Array(1..6, 1..6, {(1, 1) = -1.8174564552465413, (1, 2) = -1.5961798310015682, (1, 3) = -1.6320984766762368, (1, 4) = -1.794276191365887, (1, 5) = -1.8196239216963577, (1, 6) = -1.666791376880401, (2, 1) = 4.30366785167187, (2, 2) = 3.7285529660791323, (2, 3) = 3.8220571587930636, (2, 4) = 4.2436470774636215, (2, 5) = 4.309254368325473, (2, 6) = 3.9125540306811586, (3, 1) = 13.960107782236397, (3, 2) = 12.107313789961491, (3, 3) = 12.408495351207202, (3, 4) = 13.766676071871467, (3, 5) = 13.978102815543973, (3, 6) = 12.699960113697669, (4, 1) = 0.14695938777494755e-2, (4, 2) = 0.1890417736254646e-2, (4, 3) = 0.18192368255635083e-2, (4, 4) = 0.15081516037365184e-2, (4, 5) = 0.1464218771815906e-2, (4, 6) = 0.17469485986022493e-2, (5, 1) = 0.33726538182153962e-2, (5, 2) = 0.3616834102103271e-2, (5, 3) = 0.3575806202298267e-2, (5, 4) = 0.3396699888512361e-2, (5, 5) = 0.33718646416661044e-2, (5, 6) = 0.3534391454307867e-2, (6, 1) = -0.2756164352142944e-2, (6, 2) = -0.30603548985099053e-2, (6, 3) = -0.3009093417590724e-2, (6, 4) = -0.2785171309676129e-2, (6, 5) = -0.27538569589072135e-2, (6, 6) = -0.2957201051683882e-2}, datatype = float[8], order = C_order), Array(1..6, {(1) = 0, (2) = 0, (3) = 0, (4) = 0, (5) = 0, (6) = 0}, datatype = integer[8]), Array(1..6, {(1) = -.58658912992647, (2) = 1.567090440655536, (3) = 5.160614737970771, (4) = 0.9908058161734378e-3, (5) = .9968199579900292, (6) = 0.18025382964810283e-2}, datatype = float[8], order = C_order), Array(1..6, {(1) = -.7043836227107904, (2) = 1.8439931145266863, (3) = 6.059328015409783, (4) = 0.11097276577890164e-2, (5) = .997064605805955, (6) = 0.15985092092862313e-2}, datatype = float[8], order = C_order), Array(1..6, {(1) = 0.54780767477780046e-7, (2) = 0.1349491964841576e-6, (3) = 0.4366842834002682e-6, (4) = 0.6213567312583768e-10, (5) = 0.6205702618444775e-11, (6) = 0.24046496564789077e-10}, datatype = float[8], order = C_order), Array(1..6, {(1) = .2342493224, (2) = .1799302829, (3) = .5137385362, (4) = -0.1e-1, (5) = .99, (6) = 0.1e-1}, datatype = float[8], order = C_order), Array(1..6, {(1) = -1.7035037265233508, (2) = 4.008091885160755, (3) = 13.007711580584129, (4) = 0.16741346853412775e-2, (5) = 0.34931094516621056e-2, (6) = -0.290522356073048e-2}, datatype = float[8], order = C_order)]), ( 8 ) = ([Array(1..6, {(1) = .2342493224, (2) = .1799302829, (3) = .5137385362, (4) = -0.1e-1, (5) = .99, (6) = 0.1e-1}, datatype = float[8], order = C_order), Array(1..6, {(1) = .0, (2) = .0, (3) = .0, (4) = .0, (5) = .0, (6) = .0}, datatype = float[8], order = C_order), Array(1..6, {(1) = -1.8174564552465413, (2) = 4.30366785167187, (3) = 13.960107782236397, (4) = 0.14695938777494755e-2, (5) = 0.33726538182153962e-2, (6) = -0.2756164352142944e-2}, datatype = float[8], order = C_order), 0, 0]), ( 11 ) = (Array(1..6, 0..6, {(1, 1) = .9879192868180985, (1, 2) = -.567851996149892, (1, 3) = 1.52348122819718, (1, 4) = 5.018966795754265, (1, 5) = 0.966604406750015e-3, (1, 6) = .996775525472013, (2, 0) = .996775525472013, (2, 1) = 0.18404595686446383e-2, (2, 2) = 1.0083606305475792, (2, 3) = -.5998331149590153, (2, 4) = 1.597995931032739, (2, 5) = 5.260978676094943, (2, 6) = 0.10069043474277527e-2, (3, 0) = 0.10069043474277527e-2, (3, 1) = .9968503505551642, (3, 2) = 0.17767510111773622e-2, (3, 3) = 1.0288019742770595, (3, 4) = -.6332001141133419, (3, 5) = 1.676135536468881, (3, 6) = 5.514662587869421, (4, 0) = 5.514662587869421, (4, 1) = 0.10440945454514294e-2, (4, 2) = .9969234459541667, (4, 3) = 0.17152395072317103e-2, (4, 4) = 1.0492433180065404, (4, 5) = -.6680251208244752, (4, 6) = 1.7580738107890026, (5, 0) = 1.7580738107890026, (5, 1) = 5.780581974544998, (5, 2) = 0.1078321938733522e-2, (5, 3) = .996994851419178, (5, 4) = 0.16558498219658045e-2, (5, 5) = 1.0696846617360207, (5, 6) = -.7043836227107904, (6, 0) = -.7043836227107904, (6, 1) = 1.8439931145266863, (6, 2) = 6.059328015409783, (6, 3) = 0.11097276577890164e-2, (6, 4) = .997064605805955, (6, 5) = 0.15985092092862313e-2, (6, 6) = .0}, datatype = float[8], order = C_order)), ( 10 ) = ([proc (N, X, Y, YP) option `[Y[1] = P(t), Y[2] = Q(t), Y[3] = R(t), Y[4] = x(t), Y[5] = y(t), Y[6] = z(t)]`; YP[1] := -(1/2)*Y[3]-3.46410161513776*X(X)*Y[2]-(1.15470053837925*(-Y(X)^2+2)*X(X)+1.73205080756888*(-2*X(X)^2+Z(X)+1)*Y(X))*Y[1]; YP[2] := (2/3)*Y[3]-1.15470053837925*((1/3)*Y(X)+X(X))*Y[1]+(-1.15470053837925*(-Y(X)^2+2)*X(X)+3.46410161513776*(X(X)^2-(1/2)*Z(X)-2/3)*X(X))*Y[2]; YP[3] := Y[1]*Z(X)+(-1.15470053837925*(-Y(X)^2+2)*X(X)+1.73205080756888*(-2*X(X)^2+Z(X)+4/3)*Y(X))*Y[3]; YP[4] := 1.73205080756888*Y[4]^3*Y[5]+.288675134594813*(2*Y[5]^2-2)*Y[4]^2+.866025403784440*Y[5]*(Y[6]-2)*Y[4]-.577350269189627*Y[5]^2+.577350269189627; YP[5] := .288675134594813*(Y[5]-1)*(Y[5]+1)*(6*Y[4]^2+2*Y[4]*Y[5]+3*Y[6]-2); YP[6] := .577350269189627*Y[6]*(6*Y[4]^2*Y[5]+2*Y[4]*Y[5]^2+3*Y[5]*Y[6]-2*Y[4]-3*Y[5]); 0 end proc, -1, 0, 0, 0, 0, 0, 0]), ( 13 ) = (), ( 12 ) = (Array(1..57, 0..6, {(1, 1) = .0, (1, 2) = .2342493224, (1, 3) = .1799302829, (1, 4) = .5137385362, (1, 5) = -0.1e-1, (1, 6) = .99, (2, 0) = .99, (2, 1) = 0.1e-1, (2, 2) = 0.64545004936855275e-2, (2, 3) = .2300419710605719, (2, 4) = .18164420512820245, (2, 5) = .5214868078505989, (2, 6) = -0.9817010960328318e-2, (3, 0) = -0.9817010960328318e-2, (3, 1) = .9900734868383678, (3, 2) = 0.9891053571807595e-2, (3, 3) = 0.12909000987371055e-1, (3, 4) = .2258554096236755, (3, 5) = .18340101326334712, (3, 6) = .5293505334512683, (4, 0) = .5293505334512683, (4, 1) = -0.963654887829546e-2, (4, 2) = .9901464358841768, (4, 3) = 0.97832732764486e-2, (4, 4) = 0.1936350148105658e-1, (4, 5) = .2216890532504108, (4, 6) = .18520118185433687, (5, 0) = .18520118185433687, (5, 1) = .5373314488731611, (5, 2) = -0.9458582846696224e-2, (5, 3) = .990218851064149, (5, 4) = 0.9676647075231476e-2, (5, 5) = 0.2581800197474211e-1, (5, 6) = .21754231734236698, (6, 0) = .21754231734236698, (6, 1) = .18704519377941514, (6, 2) = .5454313160048884, (6, 3) = -0.9283082302677161e-2, (6, 4) = .990290736276485, (6, 5) = 0.9571163043176859e-2, (6, 6) = 0.4647439055831547e-1, (7, 0) = 0.4647439055831547e-1, (7, 1) = .20439654249141967, (7, 2) = .1932470634293907, (7, 3) = .5721739067233804, (7, 4) = -0.873764711557635e-2, (7, 5) = .9905172740251478, (7, 6) = 0.9241121332454215e-2, (8, 0) = 0.9241121332454215e-2, (8, 1) = 0.6713077914188882e-1, (8, 2) = .19142659263996561, (8, 3) = .19991976183931026, (8, 4) = .6002137548840873, (8, 5) = -0.8216189418555515e-2, (8, 6) = .9907385481671833, (9, 0) = .9907385481671833, (9, 1) = 0.8922278223403118e-2, (9, 2) = 0.8778716772546219e-1, (9, 3) = .17861331282447612, (9, 4) = .20708092651434187, (9, 5) = .6296142549713681, (9, 6) = -0.7717776635937243e-2, (10, 0) = -0.7717776635937243e-2, (10, 1) = .9909546807007885, (10, 2) = 0.8614265919503292e-2, (10, 3) = .10844355630903554, (10, 4) = .16593749530356997, (10, 5) = .21474919790468533, (10, 6) = .6604419364875742, (11, 0) = .6604419364875742, (11, 1) = -0.7241508776032443e-2, (11, 2) = .9911657908361249, (11, 3) = 0.8316727620072777e-2, (11, 4) = .12937487180515225, (11, 5) = .1532134446009373, (11, 6) = .22305695775391968, (12, 0) = .22305695775391968, (12, 1) = .6932071749548085, (12, 2) = -0.6780602110515472e-2, (12, 3) = .9913747069294436, (12, 4) = 0.802555894975736e-2, (12, 5) = .15030618730126896, (12, 6) = .14059064057633014, (13, 0) = .14059064057633014, (13, 1) = .2319271660714237, (13, 2) = .7275854144907258, (13, 3) = -0.6340674144478111e-2, (13, 4) = .9915787038437571, (13, 5) = 0.77444416365733204e-2, (13, 6) = .17123750279738567, (14, 0) = .17123750279738567, (14, 1) = .12804882664813438, (14, 2) = .24138260784140847, (14, 3) = .7636564614547035, (14, 4) = -0.5920887440605227e-2, (14, 5) = .991777897073362, (14, 6) = 0.7473038283094894e-2, (15, 0) = 0.7473038283094894e-2, (15, 1) = .19216881829350238, (15, 2) = .11556756924520604, (15, 3) = .2514473472807112, (15, 4) = .8015041315322827, (15, 5) = -0.5520434590675029e-2, (15, 6) = .991972399435382, (16, 0) = .991972399435382, (16, 1) = 0.7211021922350818e-2, (16, 2) = .2130181171265062, (16, 3) = .10317494576308055, (16, 4) = .2621035570104014, (16, 5) = .8410569673604387, (16, 6) = -0.51399991512683305e-2, (17, 0) = -0.51399991512683305e-2, (17, 1) = .9921615857018771, (17, 2) = 0.6959049805608757e-2, (17, 3) = .23386741595951, (17, 4) = 0.9080137433142156e-1, (17, 5) = .2734158608744957, (17, 6) = .8825507496844689, (18, 0) = .8825507496844689, (18, 1) = -0.47772357315840714e-2, (18, 2) = .9923463333316332, (18, 3) = 0.6715774612386024e-2, (18, 4) = .2547167147925138, (18, 5) = 0.784260095662413e-1, (18, 6) = .2854121087863057, (19, 0) = .2854121087863057, (19, 1) = .9260810640049988, (19, 2) = -0.4431432151446765e-2, (19, 3) = .9925267461060758, (19, 4) = 0.6480903335651844e-2, (19, 5) = .2755660136255176, (19, 6) = 0.6602769484181192e-1, (20, 0) = 0.6602769484181192e-1, (20, 1) = .29812167175990834, (20, 2) = .9717482847854149, (20, 3) = -0.4101901989273883e-2, (20, 4) = .9927029254103397, (20, 5) = 0.62541521480065625e-2, (20, 6) = .2963377210486575, (21, 0) = .2963377210486575, (21, 1) = 0.5363136453463029e-1, (21, 2) = .3115239844512747, (21, 3) = 1.0194750850357375, (21, 4) = -0.37891244661550066e-2, (21, 5) = .9928743375224754, (21, 6) = 0.6036046803078757e-2, (22, 0) = 0.6036046803078757e-2, (22, 1) = .31710942847179746, (22, 2) = 0.41169361980538305e-1, (22, 3) = .32569715100902347, (22, 4) = 1.069536411370093, (22, 5) = -0.34912151996568394e-2, (22, 6) = .9930417414423722, (23, 0) = .9930417414423722, (23, 1) = 0.5825465996974255e-2, (23, 2) = .3378811358949374, (23, 3) = 0.28619694808504743e-1, (23, 4) = .34067505054658154, (23, 5) = 1.122046757421991, (23, 6) = -0.32075689464058275e-2, (24, 0) = -0.32075689464058275e-2, (24, 1) = .9932052305205858, (24, 2) = 0.5622155526420657e-2, (24, 3) = .35865284331807734, (24, 4) = 0.1595990960349384e-1, (24, 5) = .3564933699157787, (24, 6) = 1.1771263360158977, (25, 0) = 1.1771263360158977, (25, 1) = -0.29376025300628444e-2, (25, 2) = .9933648959619354, (25, 3) = 0.5425869236039755e-2, (25, 4) = .3780899227409907, (25, 5) = 0.3993497479698814e-2, (25, 6) = .3720897166079034, (26, 0) = .3720897166079034, (26, 1) = 1.2311054664473735, (26, 2) = -0.2696873941858087e-2, (26, 3) = .9935109185048253, (26, 4) = 0.5248345298056619e-2, (26, 5) = .397527002163904, (26, 6) = -0.8108584479492708e-2, (27, 0) = -0.8108584479492708e-2, (27, 1) = .3884870439844688, (27, 2) = 1.2875525914001194, (27, 3) = -0.246718956456596e-2, (27, 4) = .9936537422859509, (27, 5) = 0.50765745538758076e-2, (27, 6) = .4169640815868174, (28, 0) = .4169640815868174, (28, 1) = -0.2036596469323423e-1, (28, 2) = .40571919462143247, (28, 3) = 1.3465807836913668, (28, 4) = -0.22481222212627616e-2, (28, 5) = .9937934370612704, (28, 6) = 0.4910373943516776e-2, (29, 0) = 0.4910373943516776e-2, (29, 1) = .43640116100973075, (29, 2) = -0.3279875611392075e-1, (29, 3) = .42382167125081455, (29, 4) = 1.4083083913938421, (29, 5) = -0.20392595109259455e-2, (29, 6) = .9939300710854048, (30, 0) = .9939300710854048, (30, 1) = 0.47495659073825735e-2, (30, 2) = .45491143444811805, (30, 3) = -0.44820635823940594e-1, (30, 4) = .4419040505816245, (30, 5) = 1.4697152115226848, (30, 6) = -0.18494782132343181e-2, (31, 0) = -0.18494782132343181e-2, (31, 1) = .9940574058527069, (31, 2) = 0.4601281096703669e-2, (31, 3) = .47342170788650534, (31, 4) = -0.5703851503536518e-1, (31, 5) = .4608430305917598, (31, 6) = 1.5337937906677732, (32, 0) = 1.5337937906677732, (32, 1) = -0.16682596095454498e-2, (32, 2) = .994182081759383, (32, 3) = 0.4457586749675641e-2, (32, 4) = .49193198132489263, (32, 5) = -0.6947111098190197e-1, (32, 6) = .4806737482572798, (33, 0) = .4806737482572798, (33, 1) = 1.6006605735100923, (33, 2) = -0.1495283606897034e-2, (33, 3) = .9943041540507893, (33, 4) = 0.4318343044245988e-2, (33, 5) = .5104422547632799, (33, 6) = -0.8213765677504237e-1, (34, 0) = -0.8213765677504237e-1, (34, 1) = .5014329607435972, (34, 2) = 1.6704371713629431, (34, 3) = -0.13302407832413619e-2, (34, 4) = .9944236768398765, (34, 5) = 0.41834142036754115e-2, (34, 6) = .5293620427048102, (35, 0) = .5293620427048102, (35, 1) = -0.9534677951685261e-1, (35, 2) = .5236509705292217, (35, 3) = 1.7448966852572163, (35, 4) = -0.1169433938920059e-2, (35, 5) = .9945432643315404, (35, 6) = 0.4049822210875571e-2, (36, 0) = 0.4049822210875571e-2, (36, 1) = .5482818306463406, (36, 2) = -.10884265049203679, (36, 3) = .546921970421614, (36, 4) = 1.8226700233244277, (36, 5) = -0.10162952917302185e-2, (36, 6) = .9946602989886729, (37, 0) = .9946602989886729, (37, 1) = 0.392046412701727e-2, (37, 2) = .567201618587871, (37, 3) = -.12264762641220278, (37, 4) = .5712908007111066, (37, 5) = 1.9039048868070536, (37, 6) = -0.8705277641523775e-3, (38, 0) = -0.8705277641523775e-3, (38, 1) = .9947748350186433, (38, 2) = 0.37952076860058325e-2, (38, 3) = .5861214065294014, (38, 4) = -.13678477408914852, (38, 5) = .5968043980443378, (38, 6) = 1.9887556822334576, (39, 0) = 1.9887556822334576, (39, 1) = -0.7318444897168429e-3, (39, 2) = .9948869254949941, (39, 3) = 0.3673924560267233e-2, (39, 4) = .6056342805983015, (39, 5) = -.15173822969663803, (39, 6) = .6243688981494192, (40, 0) = .6243688981494192, (40, 1) = 2.0802247890867167, (40, 2) = -0.5959418434413506e-3, (40, 3) = .9950000229434407, (40, 4) = 0.355286995842563e-2, (40, 5) = .6251471546672016, (40, 6) = -.16709725003017134, (41, 0) = -.16709725003017134, (41, 1) = .653259527918935, (41, 2) = 2.1758962632393546, (41, 3) = -0.4669878991351284e-3, (41, 4) = .9951106300946854, (41, 5) = 0.3435777237637031e-2, (41, 6) = .6446600287361018, (42, 0) = .6446600287361018, (42, 1) = -.18288971416024324, (42, 2) = .6835352565183088, (42, 3) = 2.275963430511651, (42, 4) = -0.3447008310403706e-3, (42, 5) = .995218801477557, (42, 6) = 0.3322518380569725e-2, (43, 0) = 0.3322518380569725e-2, (43, 1) = .664172902805002, (43, 2) = -.19914450315213802, (43, 3) = .7152578797595611, (43, 4) = 2.380628680376913, (43, 5) = -0.22880888917630805e-3, (43, 6) = .995324590446651, (44, 0) = .995324590446651, (44, 1) = 0.32129693237297856e-2, (44, 2) = .6840932625513285, (44, 3) = -.21624671715855462, (44, 4) = .7492026953426051, (44, 5) = 2.4924426415116363, (44, 6) = -0.1168217693613717e-3, (45, 0) = -0.1168217693613717e-3, (45, 1) = .9954301852399416, (45, 2) = 0.310483460561034e-2, (45, 3) = .704013622297655, (45, 4) = -.23389488251507493, (45, 5) = .7847948569407931, (45, 6) = 2.6095046231066137, (46, 0) = 2.6095046231066137, (46, 1) = -0.10959081159759578e-4, (46, 2) = .9955334057380684, (46, 3) = 0.30003174501073986e-2, (46, 4) = .7239339820439815, (46, 5) = -.252123225048779, (46, 6) = .8221098622917561, (47, 0) = .8221098622917561, (47, 1) = 2.7320612405614773, (47, 2) = 0.8903705137091036e-4, (47, 3) = .9956343050190749, (47, 4) = 0.2899298199376844e-2, (47, 5) = .7438543417903081, (47, 6) = -.27096731864504614, (48, 0) = -.27096731864504614, (48, 1) = .8612268860626413, (48, 2) = 2.8603709233982384, (48, 3) = 0.18341500433437686e-3, (48, 4) = .9957329349956677, (48, 5) = 0.2801660988365962e-2, (48, 6) = .7640460720743559, (49, 0) = .7640460720743559, (49, 1) = -.29073439837569187, (49, 2) = .9028008249282765, (49, 3) = 2.9965770264786404, (49, 4) = 0.2735900973182273e-3, (49, 5) = .9958306447281304, (49, 6) = 0.27060302381619462e-2, (50, 0) = 0.27060302381619462e-2, (50, 1) = .7842378023584038, (50, 2) = -.3112120473062878, (50, 3) = .9464017913347063, (50, 4) = 3.1392665053790916, (50, 5) = 0.3584774199777781e-3, (50, 6) = .9959261266976464, (51, 0) = .9959261266976464, (51, 1) = 0.26136463326449475e-2, (51, 2) = .8044295326424517, (51, 3) = -.3324418939267032, (51, 4) = .9921246256964269, (51, 5) = 3.2887483112560285, (51, 6) = 0.4383069259582014e-3, (52, 0) = 0.4383069259582014e-3, (52, 1) = .9960194313975269, (52, 2) = 0.2524400138847321e-2, (52, 3) = .8246212629264995, (52, 4) = -.3544673192622377, (52, 5) = 1.0400688331171313, (52, 6) = 3.445346406808373, (53, 0) = 3.445346406808373, (53, 1) = 0.5132998959894075e-3, (53, 2) = .9961106081983766, (53, 3) = 0.24381860471132406e-2, (53, 4) = .8449880119776962, (53, 5) = -.3775355071715037, (53, 6) = 1.0907848902438042, (54, 0) = 1.0907848902438042, (54, 1) = 3.6108556792520505, (54, 2) = 0.5842593671096287e-3, (54, 3) = .9962004686836631, (54, 4) = 0.23541925639434195e-2, (54, 5) = .865354761028893, (54, 6) = -.4015073284744735, (55, 0) = -.4015073284744735, (55, 1) = 1.143979084502326, (55, 2) = 3.784314630295771, (55, 3) = 0.6507239165973375e-3, (55, 4) = .9962882616892705, (55, 5) = 0.22730788940944855e-2, (55, 6) = .8857215100800896, (56, 0) = .8857215100800896, (56, 1) = -.4264331477087061, (56, 2) = 1.199768991342326, (56, 3) = 3.9661054923906045, (56, 4) = 0.7128948846262493e-3, (56, 5) = .9963740345004575, (56, 6) = 0.21947471776414245e-2, (57, 0) = 0.21947471776414245e-2, (57, 1) = .9060882591312863, (57, 2) = -.4523655540495244, (57, 3) = 1.2582780026726232, (57, 4) = 4.156629239071036, (57, 5) = 0.7709658694348732e-3, (57, 6) = .9964578333425145}, datatype = float[8], order = C_order)), ( 15 ) = ("rkf45"), ( 14 ) = ([0, 0]), ( 19 ) = (0), ( 16 ) = ([0, 0, 0, []]), ( 17 ) = ([proc (N, X, Y, YP) option `[Y[1] = P(t), Y[2] = Q(t), Y[3] = R(t), Y[4] = x(t), Y[5] = y(t), Y[6] = z(t)]`; YP[1] := -(1/2)*Y[3]-3.46410161513776*X(X)*Y[2]-(1.15470053837925*(-Y(X)^2+2)*X(X)+1.73205080756888*(-2*X(X)^2+Z(X)+1)*Y(X))*Y[1]; YP[2] := (2/3)*Y[3]-1.15470053837925*((1/3)*Y(X)+X(X))*Y[1]+(-1.15470053837925*(-Y(X)^2+2)*X(X)+3.46410161513776*(X(X)^2-(1/2)*Z(X)-2/3)*X(X))*Y[2]; YP[3] := Y[1]*Z(X)+(-1.15470053837925*(-Y(X)^2+2)*X(X)+1.73205080756888*(-2*X(X)^2+Z(X)+4/3)*Y(X))*Y[3]; YP[4] := 1.73205080756888*Y[4]^3*Y[5]+.288675134594813*(2*Y[5]^2-2)*Y[4]^2+.866025403784440*Y[5]*(Y[6]-2)*Y[4]-.577350269189627*Y[5]^2+.577350269189627; YP[5] := .288675134594813*(Y[5]-1)*(Y[5]+1)*(6*Y[4]^2+2*Y[4]*Y[5]+3*Y[6]-2); YP[6] := .577350269189627*Y[6]*(6*Y[4]^2*Y[5]+2*Y[4]*Y[5]^2+3*Y[5]*Y[6]-2*Y[4]-3*Y[5]); 0 end proc, -1, 0, 0, 0, 0, 0, 0]), ( 22 ) = (0), ( 23 ) = (0), ( 20 ) = ([]), ( 21 ) = (0), ( 24 ) = (0) ] )), ( 4 ) = (3) ] ); _y0 := Array(0..6, {(1) = 0., (2) = .2342493224, (3) = .1799302829, (4) = .5137385362, (5) = -0.1e-1, (6) = .99}); _vmap := array( 1 .. 6, [( 1 ) = (1), ( 2 ) = (2), ( 3 ) = (3), ( 4 ) = (4), ( 5 ) = (5), ( 6 ) = (6) ] ); _x0 := _dtbl[1][5][5]; _n := _dtbl[1][4][1]; _ne := _dtbl[1][4][3]; _nd := _dtbl[1][4][4]; _nv := _dtbl[1][4][16]; if not type(_xout, 'numeric') then if member(_xout, ["start", "left", "right"]) then if _Env_smart_dsolve_numeric = true or _dtbl[1][4][10] = 1 then if _xout = "left" then if type(_dtbl[2], 'table') then return _dtbl[2][5][1] end if elif _xout = "right" then if type(_dtbl[3], 'table') then return _dtbl[3][5][1] end if end if end if; return _dtbl[1][5][5] elif _xout = "method" then return _dtbl[1][15] elif _xout = "storage" then return evalb(_dtbl[1][4][10] = 1) elif _xout = "leftdata" then if not type(_dtbl[2], 'array') then return NULL else return eval(_dtbl[2]) end if elif _xout = "rightdata" then if not type(_dtbl[3], 'array') then return NULL else return eval(_dtbl[3]) end if elif _xout = "enginedata" then return eval(_dtbl[1]) elif _xout = "enginereset" then _dtbl[2] := evaln(_dtbl[2]); _dtbl[3] := evaln(_dtbl[3]); return NULL elif _xout = "initial" then return procname(_y0[0]) elif _xout = "laxtol" then return _dtbl[`if`(member(_dtbl[4], {2, 3}), _dtbl[4], 1)][5][18] elif _xout = "numfun" then return `if`(member(_dtbl[4], {2, 3}), _dtbl[_dtbl[4]][4][18], 0) elif _xout = "parameters" then return [seq(_y0[_n+_i], _i = 1 .. nops(_pars))] elif _xout = "initial_and_parameters" then return procname(_y0[0]), [seq(_y0[_n+_i], _i = 1 .. nops(_pars))] elif _xout = "last" then if _dtbl[4] <> 2 and _dtbl[4] <> 3 or _x0-_dtbl[_dtbl[4]][5][1] = 0. then error "no information is available on last computed point" else _xout := _dtbl[_dtbl[4]][5][1] end if elif _xout = "function" then if _dtbl[1][4][33]-2. = 0 then return eval(_dtbl[1][10], 1) else return eval(_dtbl[1][10][1], 1) end if elif _xout = "map" then return copy(_vmap) elif type(_xin, `=`) and type(rhs(_xin), 'list') and member(lhs(_xin), {"initial", "parameters", "initial_and_parameters"}) then _ini, _par := [], []; if lhs(_xin) = "initial" then _ini := rhs(_xin) elif lhs(_xin) = "parameters" then _par := rhs(_xin) elif select(type, rhs(_xin), `=`) <> [] then _par, _ini := selectremove(type, rhs(_xin), `=`) elif nops(rhs(_xin)) < nops(_pars)+1 then error "insufficient data for specification of initial and parameters" else _par := rhs(_xin)[-nops(_pars) .. -1]; _ini := rhs(_xin)[1 .. -nops(_pars)-1] end if; _xout := lhs(_xout); if _par <> [] then `dsolve/numeric/process_parameters`(_n, _pars, _par, _y0) end if; if _ini <> [] then `dsolve/numeric/process_initial`(_n-_ne, _ini, _y0, _pars, _vmap) end if; `dsolve/numeric/SC/reinitialize`(_dtbl, _y0, _n, procname, _pars); if _Env_smart_dsolve_numeric = true and type(_y0[0], 'numeric') and _dtbl[1][4][10] <> 1 then procname("right") := _y0[0]; procname("left") := _y0[0] end if; if _xout = "initial" then return [_y0[0], seq(_y0[_vmap[_i]], _i = 1 .. _n-_ne)] elif _xout = "parameters" then return [seq(_y0[_n+_i], _i = 1 .. nops(_pars))] else return [_y0[0], seq(_y0[_vmap[_i]], _i = 1 .. _n-_ne)], [seq(_y0[_n+_i], _i = 1 .. nops(_pars))] end if elif _xin = "eventstop" then if _nv = 0 then error "this solution has no events" end if; _i := _dtbl[4]; if _i <> 2 and _i <> 3 then return 0 end if; if _dtbl[_i][4][10] = 1 and assigned(_dtbl[5-_i]) and _dtbl[_i][4][9] < 100 and 100 <= _dtbl[5-_i][4][9] then _i := 5-_i; _dtbl[4] := _i; _j := round(_dtbl[_i][4][17]); return round(_dtbl[_i][3][1][_j, 1]) elif 100 <= _dtbl[_i][4][9] then _j := round(_dtbl[_i][4][17]); return round(_dtbl[_i][3][1][_j, 1]) else return 0 end if elif _xin = "eventstatus" then if _nv = 0 then error "this solution has no events" end if; _i := [selectremove(proc (a) options operator, arrow; _dtbl[1][3][1][a, 7] = 1 end proc, {seq(_j, _j = 1 .. round(_dtbl[1][3][1][_nv+1, 1]))})]; return ':-enabled' = _i[1], ':-disabled' = _i[2] elif _xin = "eventclear" then if _nv = 0 then error "this solution has no events" end if; _i := _dtbl[4]; if _i <> 2 and _i <> 3 then error "no events to clear" end if; if _dtbl[_i][4][10] = 1 and assigned(_dtbl[5-_i]) and _dtbl[_i][4][9] < 100 and 100 < _dtbl[5-_i][4][9] then _dtbl[4] := 5-_i; _i := 5-_i end if; if _dtbl[_i][4][9] < 100 then error "no events to clear" elif _nv < _dtbl[_i][4][9]-100 then error "event error condition cannot be cleared" else _j := _dtbl[_i][4][9]-100; if irem(round(_dtbl[_i][3][1][_j, 4]), 2) = 1 then error "retriggerable events cannot be cleared" end if; _j := round(_dtbl[_i][3][1][_j, 1]); for _k to _nv do if _dtbl[_i][3][1][_k, 1] = _j then if _dtbl[_i][3][1][_k, 2] = 3 then error "range events cannot be cleared" end if; _dtbl[_i][3][1][_k, 8] := _dtbl[_i][3][1][_nv+1, 8] end if end do; _dtbl[_i][4][17] := 0; _dtbl[_i][4][9] := 0; if _dtbl[1][4][10] = 1 then if _i = 2 then try procname(procname("left")) catch: end try else try procname(procname("right")) catch: end try end if end if end if; return elif type(_xin, `=`) and member(lhs(_xin), {"eventdisable", "eventenable"}) then if _nv = 0 then error "this solution has no events" end if; if type(rhs(_xin), {('list')('posint'), ('set')('posint')}) then _i := {op(rhs(_xin))} elif type(rhs(_xin), 'posint') then _i := {rhs(_xin)} else error "event identifiers must be integers in the range 1..%1", round(_dtbl[1][3][1][_nv+1, 1]) end if; if select(proc (a) options operator, arrow; _nv < a end proc, _i) <> {} then error "event identifiers must be integers in the range 1..%1", round(_dtbl[1][3][1][_nv+1, 1]) end if; _k := {}; for _j to _nv do if member(round(_dtbl[1][3][1][_j, 1]), _i) then _k := `union`(_k, {_j}) end if end do; _i := _k; if lhs(_xin) = "eventdisable" then _dtbl[4] := 0; _j := [evalb(assigned(_dtbl[2]) and member(_dtbl[2][4][17], _i)), evalb(assigned(_dtbl[3]) and member(_dtbl[3][4][17], _i))]; for _k in _i do _dtbl[1][3][1][_k, 7] := 0; if assigned(_dtbl[2]) then _dtbl[2][3][1][_k, 7] := 0 end if; if assigned(_dtbl[3]) then _dtbl[3][3][1][_k, 7] := 0 end if end do; if _j[1] then for _k to _nv+1 do if _k <= _nv and not type(_dtbl[2][3][4][_k, 1], 'undefined') then userinfo(3, {'events', 'eventreset'}, `reinit #2, event code `, _k, ` to defined init `, _dtbl[2][3][4][_k, 1]); _dtbl[2][3][1][_k, 8] := _dtbl[2][3][4][_k, 1] elif _dtbl[2][3][1][_k, 2] = 0 and irem(iquo(round(_dtbl[2][3][1][_k, 4]), 32), 2) = 1 then userinfo(3, {'events', 'eventreset'}, `reinit #2, event code `, _k, ` to rate hysteresis init `, _dtbl[2][5][24]); _dtbl[2][3][1][_k, 8] := _dtbl[2][5][24] elif _dtbl[2][3][1][_k, 2] = 0 and irem(iquo(round(_dtbl[2][3][1][_k, 4]), 2), 2) = 0 then userinfo(3, {'events', 'eventreset'}, `reinit #2, event code `, _k, ` to initial init `, _x0); _dtbl[2][3][1][_k, 8] := _x0 else userinfo(3, {'events', 'eventreset'}, `reinit #2, event code `, _k, ` to fireinitial init `, _x0-1); _dtbl[2][3][1][_k, 8] := _x0-1 end if end do; _dtbl[2][4][17] := 0; _dtbl[2][4][9] := 0; if _dtbl[1][4][10] = 1 then procname(procname("left")) end if end if; if _j[2] then for _k to _nv+1 do if _k <= _nv and not type(_dtbl[3][3][4][_k, 2], 'undefined') then userinfo(3, {'events', 'eventreset'}, `reinit #3, event code `, _k, ` to defined init `, _dtbl[3][3][4][_k, 2]); _dtbl[3][3][1][_k, 8] := _dtbl[3][3][4][_k, 2] elif _dtbl[3][3][1][_k, 2] = 0 and irem(iquo(round(_dtbl[3][3][1][_k, 4]), 32), 2) = 1 then userinfo(3, {'events', 'eventreset'}, `reinit #3, event code `, _k, ` to rate hysteresis init `, _dtbl[3][5][24]); _dtbl[3][3][1][_k, 8] := _dtbl[3][5][24] elif _dtbl[3][3][1][_k, 2] = 0 and irem(iquo(round(_dtbl[3][3][1][_k, 4]), 2), 2) = 0 then userinfo(3, {'events', 'eventreset'}, `reinit #3, event code `, _k, ` to initial init `, _x0); _dtbl[3][3][1][_k, 8] := _x0 else userinfo(3, {'events', 'eventreset'}, `reinit #3, event code `, _k, ` to fireinitial init `, _x0+1); _dtbl[3][3][1][_k, 8] := _x0+1 end if end do; _dtbl[3][4][17] := 0; _dtbl[3][4][9] := 0; if _dtbl[1][4][10] = 1 then procname(procname("right")) end if end if else for _k in _i do _dtbl[1][3][1][_k, 7] := 1 end do; _dtbl[2] := evaln(_dtbl[2]); _dtbl[3] := evaln(_dtbl[3]); _dtbl[4] := 0; if _dtbl[1][4][10] = 1 then if _x0 <= procname("right") then try procname(procname("right")) catch: end try end if; if procname("left") <= _x0 then try procname(procname("left")) catch: end try end if end if end if; return elif type(_xin, `=`) and lhs(_xin) = "eventfired" then if not type(rhs(_xin), 'list') then error "'eventfired' must be specified as a list" end if; if _nv = 0 then error "this solution has no events" end if; if _dtbl[4] <> 2 and _dtbl[4] <> 3 then error "'direction' must be set prior to calling/setting 'eventfired'" end if; _i := _dtbl[4]; _val := NULL; if not assigned(_EnvEventRetriggerWarned) then _EnvEventRetriggerWarned := false end if; for _k in rhs(_xin) do if type(_k, 'integer') then _src := _k elif type(_k, 'integer' = 'anything') and type(evalf(rhs(_k)), 'numeric') then _k := lhs(_k) = evalf[max(Digits, 18)](rhs(_k)); _src := lhs(_k) else error "'eventfired' entry is not valid: %1", _k end if; if _src < 1 or round(_dtbl[1][3][1][_nv+1, 1]) < _src then error "event identifiers must be integers in the range 1..%1", round(_dtbl[1][3][1][_nv+1, 1]) end if; _src := {seq(`if`(_dtbl[1][3][1][_j, 1]-_src = 0., _j, NULL), _j = 1 .. _nv)}; if nops(_src) <> 1 then error "'eventfired' can only be set/queried for root-finding events and time/interval events" end if; _src := _src[1]; if _dtbl[1][3][1][_src, 2] <> 0. and _dtbl[1][3][1][_src, 2]-2. <> 0. then error "'eventfired' can only be set/queried for root-finding events and time/interval events" elif irem(round(_dtbl[1][3][1][_src, 4]), 2) = 1 then if _EnvEventRetriggerWarned = false then WARNING(`'eventfired' has no effect on events that retrigger`) end if; _EnvEventRetriggerWarned := true end if; if _dtbl[_i][3][1][_src, 2] = 0 and irem(iquo(round(_dtbl[_i][3][1][_src, 4]), 32), 2) = 1 then _val := _val, undefined elif type(_dtbl[_i][3][4][_src, _i-1], 'undefined') or _i = 2 and _dtbl[2][3][1][_src, 8] < _dtbl[2][3][4][_src, 1] or _i = 3 and _dtbl[3][3][4][_src, 2] < _dtbl[3][3][1][_src, 8] then _val := _val, _dtbl[_i][3][1][_src, 8] else _val := _val, _dtbl[_i][3][4][_src, _i-1] end if; if type(_k, `=`) then if _dtbl[_i][3][1][_src, 2] = 0 and irem(iquo(round(_dtbl[_i][3][1][_src, 4]), 32), 2) = 1 then error "cannot set event code for a rate hysteresis event" end if; userinfo(3, {'events', 'eventreset'}, `manual set event code `, _src, ` to value `, rhs(_k)); _dtbl[_i][3][1][_src, 8] := rhs(_k); _dtbl[_i][3][4][_src, _i-1] := rhs(_k) end if end do; return [_val] elif type(_xin, `=`) and lhs(_xin) = "direction" then if not member(rhs(_xin), {-1, 1, ':-left', ':-right'}) then error "'direction' must be specified as either '1' or 'right' (positive) or '-1' or 'left' (negative)" end if; _src := `if`(_dtbl[4] = 2, -1, `if`(_dtbl[4] = 3, 1, undefined)); _i := `if`(member(rhs(_xin), {1, ':-right'}), 3, 2); _dtbl[4] := _i; _dtbl[_i] := `dsolve/numeric/SC/IVPdcopy`(_dtbl[1], `if`(assigned(_dtbl[_i]), _dtbl[_i], NULL)); if 0 < _nv then for _j to _nv+1 do if _j <= _nv and not type(_dtbl[_i][3][4][_j, _i-1], 'undefined') then userinfo(3, {'events', 'eventreset'}, `reinit #4, event code `, _j, ` to defined init `, _dtbl[_i][3][4][_j, _i-1]); _dtbl[_i][3][1][_j, 8] := _dtbl[_i][3][4][_j, _i-1] elif _dtbl[_i][3][1][_j, 2] = 0 and irem(iquo(round(_dtbl[_i][3][1][_j, 4]), 32), 2) = 1 then userinfo(3, {'events', 'eventreset'}, `reinit #4, event code `, _j, ` to rate hysteresis init `, _dtbl[_i][5][24]); _dtbl[_i][3][1][_j, 8] := _dtbl[_i][5][24] elif _dtbl[_i][3][1][_j, 2] = 0 and irem(iquo(round(_dtbl[_i][3][1][_j, 4]), 2), 2) = 0 then userinfo(3, {'events', 'eventreset'}, `reinit #4, event code `, _j, ` to initial init `, _x0); _dtbl[_i][3][1][_j, 8] := _x0 else userinfo(3, {'events', 'eventreset'}, `reinit #4, event code `, _j, ` to fireinitial init `, _x0-2*_i+5.0); _dtbl[_i][3][1][_j, 8] := _x0-2*_i+5.0 end if end do end if; return _src elif _xin = "eventcount" then if _dtbl[1][3][1] = 0 or _dtbl[4] <> 2 and _dtbl[4] <> 3 then return 0 else return round(_dtbl[_dtbl[4]][3][1][_nv+1, 12]) end if else return "procname" end if end if; if _xout = _x0 then return [_x0, seq(evalf(_dtbl[1][6][_vmap[_i]]), _i = 1 .. _n-_ne)] end if; _i := `if`(_x0 <= _xout, 3, 2); if _xin = "last" and 0 < _dtbl[_i][4][9] and _dtbl[_i][4][9] < 100 then _dat := eval(_dtbl[_i], 2); _j := _dat[4][20]; return [_dat[11][_j, 0], seq(_dat[11][_j, _vmap[_i]], _i = 1 .. _n-_ne-_nd), seq(_dat[8][1][_vmap[_i]], _i = _n-_ne-_nd+1 .. _n-_ne)] end if; if not type(_dtbl[_i], 'array') then _dtbl[_i] := `dsolve/numeric/SC/IVPdcopy`(_dtbl[1], `if`(assigned(_dtbl[_i]), _dtbl[_i], NULL)); if 0 < _nv then for _j to _nv+1 do if _j <= _nv and not type(_dtbl[_i][3][4][_j, _i-1], 'undefined') then userinfo(3, {'events', 'eventreset'}, `reinit #5, event code `, _j, ` to defined init `, _dtbl[_i][3][4][_j, _i-1]); _dtbl[_i][3][1][_j, 8] := _dtbl[_i][3][4][_j, _i-1] elif _dtbl[_i][3][1][_j, 2] = 0 and irem(iquo(round(_dtbl[_i][3][1][_j, 4]), 32), 2) = 1 then userinfo(3, {'events', 'eventreset'}, `reinit #5, event code `, _j, ` to rate hysteresis init `, _dtbl[_i][5][24]); _dtbl[_i][3][1][_j, 8] := _dtbl[_i][5][24] elif _dtbl[_i][3][1][_j, 2] = 0 and irem(iquo(round(_dtbl[_i][3][1][_j, 4]), 2), 2) = 0 then userinfo(3, {'events', 'eventreset'}, `reinit #5, event code `, _j, ` to initial init `, _x0); _dtbl[_i][3][1][_j, 8] := _x0 else userinfo(3, {'events', 'eventreset'}, `reinit #5, event code `, _j, ` to fireinitial init `, _x0-2*_i+5.0); _dtbl[_i][3][1][_j, 8] := _x0-2*_i+5.0 end if end do end if end if; if _xin <> "last" then if 0 < 0 then if `dsolve/numeric/checkglobals`(op(_dtbl[1][14]), _pars, _n, _y0) then `dsolve/numeric/SC/reinitialize`(_dtbl, _y0, _n, procname, _pars, _i) end if end if; if _dtbl[1][4][7] = 0 then error "parameters must be initialized before solution can be computed" end if end if; _dat := eval(_dtbl[_i], 2); _dtbl[4] := _i; try _src := `dsolve/numeric/SC/IVPrun`(_dat, _xout) catch: userinfo(2, `dsolve/debug`, print(`Exception in solnproc:`, [lastexception][2 .. -1])); error end try; if _src = 0 and 100 < _dat[4][9] then _val := _dat[3][1][_nv+1, 8] else _val := _dat[11][_dat[4][20], 0] end if; if _src <> 0 or _dat[4][9] <= 0 then _dtbl[1][5][1] := _xout else _dtbl[1][5][1] := _val end if; if _i = 3 and _val < _xout then Rounding := -infinity; if _dat[4][9] = 1 then error "cannot evaluate the solution further right of %1, probably a singularity", evalf[8](_val) elif _dat[4][9] = 2 then error "cannot evaluate the solution further right of %1, maxfun limit exceeded (see ?dsolve,maxfun for details)", evalf[8](_val) elif _dat[4][9] = 3 then if _dat[4][25] = 3 then error "cannot evaluate the solution past the initial point, problem may be initially singular or improperly set up" else error "cannot evaluate the solution past the initial point, problem may be complex, initially singular or improperly set up" end if elif _dat[4][9] = 4 then error "cannot evaluate the solution further right of %1, accuracy goal cannot be achieved with specified 'minstep'", evalf[8](_val) elif _dat[4][9] = 5 then error "cannot evaluate the solution further right of %1, too many step failures, tolerances may be too loose for problem", evalf[8](_val) elif _dat[4][9] = 6 then error "cannot evaluate the solution further right of %1, cannot downgrade delay storage for problems with delay derivative order > 1, try increasing delaypts", evalf[8](_val) elif _dat[4][9] = 10 then error "cannot evaluate the solution further right of %1, interrupt requested", evalf[8](_val) elif 100 < _dat[4][9] then if _dat[4][9]-100 = _nv+1 then error "constraint projection failure on event at t=%1", evalf[8](_val) elif _dat[4][9]-100 = _nv+2 then error "index-1 and derivative evaluation failure on event at t=%1", evalf[8](_val) elif _dat[4][9]-100 = _nv+3 then error "maximum number of event iterations reached (%1) at t=%2", round(_dat[3][1][_nv+1, 3]), evalf[8](_val) else if _Env_dsolve_nowarnstop <> true then `dsolve/numeric/warning`(StringTools:-FormatMessage("cannot evaluate the solution further right of %1, event #%2 triggered a halt", evalf[8](_val), round(_dat[3][1][_dat[4][9]-100, 1]))) end if; Rounding := 'nearest'; _xout := _val end if else error "cannot evaluate the solution further right of %1", evalf[8](_val) end if elif _i = 2 and _xout < _val then Rounding := infinity; if _dat[4][9] = 1 then error "cannot evaluate the solution further left of %1, probably a singularity", evalf[8](_val) elif _dat[4][9] = 2 then error "cannot evaluate the solution further left of %1, maxfun limit exceeded (see ?dsolve,maxfun for details)", evalf[8](_val) elif _dat[4][9] = 3 then if _dat[4][25] = 3 then error "cannot evaluate the solution past the initial point, problem may be initially singular or improperly set up" else error "cannot evaluate the solution past the initial point, problem may be complex, initially singular or improperly set up" end if elif _dat[4][9] = 4 then error "cannot evaluate the solution further left of %1, accuracy goal cannot be achieved with specified 'minstep'", evalf[8](_val) elif _dat[4][9] = 5 then error "cannot evaluate the solution further left of %1, too many step failures, tolerances may be too loose for problem", evalf[8](_val) elif _dat[4][9] = 6 then error "cannot evaluate the solution further left of %1, cannot downgrade delay storage for problems with delay derivative order > 1, try increasing delaypts", evalf[8](_val) elif _dat[4][9] = 10 then error "cannot evaluate the solution further right of %1, interrupt requested", evalf[8](_val) elif 100 < _dat[4][9] then if _dat[4][9]-100 = _nv+1 then error "constraint projection failure on event at t=%1", evalf[8](_val) elif _dat[4][9]-100 = _nv+2 then error "index-1 and derivative evaluation failure on event at t=%1", evalf[8](_val) elif _dat[4][9]-100 = _nv+3 then error "maximum number of event iterations reached (%1) at t=%2", round(_dat[3][1][_nv+1, 3]), evalf[8](_val) else if _Env_dsolve_nowarnstop <> true then `dsolve/numeric/warning`(StringTools:-FormatMessage("cannot evaluate the solution further left of %1, event #%2 triggered a halt", evalf[8](_val), round(_dat[3][1][_dat[4][9]-100, 1]))) end if; Rounding := 'nearest'; _xout := _val end if else error "cannot evaluate the solution further left of %1", evalf[8](_val) end if end if; if _EnvInFsolve = true then _dig := _dat[4][26]; _dat[4][26] := _EnvDSNumericSaveDigits; _Env_dsolve_SC_native := true; if _dat[4][25] = 1 then _i := 1; _dat[4][25] := 2 else _i := _dat[4][25] end if; _val := `dsolve/numeric/SC/IVPval`(_dat, _xout, _src); _dat[4][25] := _i; _dat[4][26] := _dig; [_xout, seq(_val[_vmap[_i]], _i = 1 .. _n-_ne)] else Digits := _dat[4][26]; _val := `dsolve/numeric/SC/IVPval`(eval(_dat, 2), _xout, _src); [_xout, seq(_val[_vmap[_i]], _i = 1 .. _n-_ne)] end if end proc, (2) = Array(0..0, {}), (3) = [t, P(t), Q(t), R(t), x(t), y(t), z(t)], (4) = []}); _vars := _dat[3]; _pars := map(rhs, _dat[4]); _n := nops(_vars)-1; _solnproc := _dat[1]; if not type(_xout, 'numeric') then if member(x_rkf45, ["start", 'start', "method", 'method', "left", 'left', "right", 'right', "leftdata", "rightdata", "enginedata", "eventstop", 'eventstop', "eventclear", 'eventclear', "eventstatus", 'eventstatus', "eventcount", 'eventcount', "laxtol", 'laxtol', "numfun", 'numfun', NULL]) then _res := _solnproc(convert(x_rkf45, 'string')); if 1 < nops([_res]) then return _res elif type(_res, 'array') then return eval(_res, 1) elif _res <> "procname" then return _res end if elif member(x_rkf45, ["last", 'last', "initial", 'initial', "parameters", 'parameters', "initial_and_parameters", 'initial_and_parameters', NULL]) then _xout := convert(x_rkf45, 'string'); _res := _solnproc(_xout); if _xout = "parameters" then return [seq(_pars[_i] = _res[_i], _i = 1 .. nops(_pars))] elif _xout = "initial_and_parameters" then return [seq(_vars[_i+1] = [_res][1][_i+1], _i = 0 .. _n), seq(_pars[_i] = [_res][2][_i], _i = 1 .. nops(_pars))] else return [seq(_vars[_i+1] = _res[_i+1], _i = 0 .. _n)] end if elif type(_xout, `=`) and member(lhs(_xout), ["initial", 'initial', "parameters", 'parameters', "initial_and_parameters", 'initial_and_parameters', NULL]) then _xout := convert(lhs(x_rkf45), 'string') = rhs(x_rkf45); if type(rhs(_xout), 'list') then _res := _solnproc(_xout) else error "initial and/or parameter values must be specified in a list" end if; if lhs(_xout) = "initial" then return [seq(_vars[_i+1] = _res[_i+1], _i = 0 .. _n)] elif lhs(_xout) = "parameters" then return [seq(_pars[_i] = _res[_i], _i = 1 .. nops(_pars))] else return [seq(_vars[_i+1] = [_res][1][_i+1], _i = 0 .. _n), seq(_pars[_i] = [_res][2][_i], _i = 1 .. nops(_pars))] end if elif type(_xout, `=`) and member(lhs(_xout), ["eventdisable", 'eventdisable', "eventenable", 'eventenable', "eventfired", 'eventfired', "direction", 'direction', NULL]) then return _solnproc(convert(lhs(x_rkf45), 'string') = rhs(x_rkf45)) elif _xout = "solnprocedure" then return eval(_solnproc) elif _xout = "sysvars" then return _vars end if; if procname <> unknown then return ('procname')(x_rkf45) else _ndsol; _ndsol := pointto(_dat[2][0]); return ('_ndsol')(x_rkf45) end if end if; try _res := _solnproc(_xout); [seq(_vars[_i+1] = _res[_i+1], _i = 0 .. _n)] catch: error end try end proc

 (6)
> 

plots:-display(
  plots:-odeplot(SOL, [t, P(t)], t=0..1, color=red, thickness=3, legend="two systems"),
  plots:-odeplot(soln, [t, P(t)], t=0..1, color=blue, thickness=7, transparency=0.8, legend="one system")
)

 
 

 

Download known_option.mw

Do you mean something like this?...

mmcdara 7284
December 16 2024
1 2

(sorry the content cannot be loaded, likely due to some plot command)

Maybe_this.mw

Illustration : epsilon_r for N=1..10 and P from 0.1 to 10

Nice (?) splitting of a latex expression...

mmcdara 7284
December 15 2024
0 4

@salim-barzani 

As an example I used your lengthy eq3 expression without any simplifications, just to work on a big expression.

The idea is to split the Maple expression on several lines by adding break lines between operators in order ot introduce a smart breaking.
Your lengthy expression is then made of several sub-expressions.
It suffices now to generate the LaTeX code for each sub-expression, to concatenate to of those code the "\\\\" breaking line command, finally to concatenate all those codes.

NOTE: the `#mo("1",mathcolor="white")` trick enables having a subexpression beginning with a '+' character (which would be automaticlly removed by Maple)

> 

restart
> 

eq3 := (4*(-(2*(lambda*B[1]^2-lambda*B[2]^2-mu^2/lambda))*lambda+4*mu^2))*alpha[1]^2*a[5]*alpha[0]-12*mu^2*alpha[1]^2*a[5]*alpha[0]+(3*(lambda*B[1]^2-lambda*B[2]^2-mu^2/lambda))*alpha[1]^2*alpha[0]*a[2]-(lambda*B[1]^2-lambda*B[2]^2-mu^2/lambda)*k^2*a[1]*alpha[1]^2+(1/2)*(lambda*B[1]^2-lambda*B[2]^2-mu^2/lambda)*alpha[1]^2*lambda*a[1]+(5*(-(2*(lambda*B[1]^2-lambda*B[2]^2-mu^2/lambda))*lambda+4*mu^2))*alpha[1]^4*alpha[0]*a[4]+(10*(lambda*B[1]^2-lambda*B[2]^2-mu^2/lambda))*alpha[1]^2*alpha[0]^3*a[4]+(6*(lambda*B[1]^2-lambda*B[2]^2-mu^2/lambda))*alpha[1]^2*alpha[0]^2*a[3]-6*lambda*beta[0]^2*alpha[1]^2*a[3]-2*lambda*beta[0]^2*a[5]*alpha[0]+6*mu*beta[0]*alpha[1]^2*a[2]+3*mu*beta[0]*a[5]*alpha[0]^2-9*mu^2*alpha[1]^2*a[1]*(1/4)+3*mu*a[1]*alpha[0]*beta[0]*(1/2)+(4*(lambda*B[1]^2-lambda*B[2]^2-mu^2/lambda))*alpha[1]^2*lambda*a[5]*alpha[0]-20*mu*beta[0]*lambda*alpha[1]^4*a[4]+24*mu*beta[0]*alpha[1]^2*alpha[0]*a[3]-30*lambda*beta[0]^2*alpha[1]^2*alpha[0]*a[4]+60*mu*beta[0]*alpha[1]^2*alpha[0]^2*a[4]-7*mu*beta[0]*lambda*a[5]*alpha[1]^2-w*beta[0]^2-(1/4)*lambda*beta[0]^2*a[1]+3*beta[0]^2*alpha[0]*a[2]-k^2*a[1]*beta[0]^2+10*beta[0]^2*alpha[0]^3*a[4]+6*beta[0]^2*alpha[0]^2*a[3]-(lambda*B[1]^2-lambda*B[2]^2-mu^2/lambda)*w*alpha[1]^2+(1/4)*(3*(-(2*(lambda*B[1]^2-lambda*B[2]^2-mu^2/lambda))*lambda+4*mu^2))*alpha[1]^2*a[1]+(-(2*(lambda*B[1]^2-lambda*B[2]^2-mu^2/lambda))*lambda+4*mu^2)*alpha[1]^4*a[3] = 0;

-(9/4)*mu^2*alpha[1]^2*a[1]+3*beta[0]^2*alpha[0]*a[2]+10*beta[0]^2*alpha[0]^3*a[4]+6*beta[0]^2*alpha[0]^2*a[3]-(1/4)*lambda*beta[0]^2*a[1]+(1/4)*(-6*(lambda*B[1]^2-lambda*B[2]^2-mu^2/lambda)*lambda+12*mu^2)*alpha[1]^2*a[1]-w*beta[0]^2+4*(lambda*B[1]^2-lambda*B[2]^2-mu^2/lambda)*alpha[1]^2*lambda*a[5]*alpha[0]-20*mu*beta[0]*lambda*alpha[1]^4*a[4]+24*mu*beta[0]*alpha[1]^2*alpha[0]*a[3]-30*lambda*beta[0]^2*alpha[1]^2*alpha[0]*a[4]+60*mu*beta[0]*alpha[1]^2*alpha[0]^2*a[4]-7*mu*beta[0]*lambda*a[5]*alpha[1]^2-(lambda*B[1]^2-lambda*B[2]^2-mu^2/lambda)*w*alpha[1]^2+(-2*(lambda*B[1]^2-lambda*B[2]^2-mu^2/lambda)*lambda+4*mu^2)*alpha[1]^4*a[3]+4*(-2*(lambda*B[1]^2-lambda*B[2]^2-mu^2/lambda)*lambda+4*mu^2)*alpha[1]^2*a[5]*alpha[0]-12*mu^2*alpha[1]^2*a[5]*alpha[0]+3*(lambda*B[1]^2-lambda*B[2]^2-mu^2/lambda)*alpha[1]^2*alpha[0]*a[2]+(1/2)*(lambda*B[1]^2-lambda*B[2]^2-mu^2/lambda)*alpha[1]^2*lambda*a[1]+5*(-2*(lambda*B[1]^2-lambda*B[2]^2-mu^2/lambda)*lambda+4*mu^2)*alpha[1]^4*alpha[0]*a[4]+10*(lambda*B[1]^2-lambda*B[2]^2-mu^2/lambda)*alpha[1]^2*alpha[0]^3*a[4]+6*(lambda*B[1]^2-lambda*B[2]^2-mu^2/lambda)*alpha[1]^2*alpha[0]^2*a[3]-6*lambda*beta[0]^2*alpha[1]^2*a[3]-2*lambda*beta[0]^2*a[5]*alpha[0]+6*mu*beta[0]*alpha[1]^2*a[2]+3*mu*beta[0]*a[5]*alpha[0]^2+(3/2)*mu*a[1]*alpha[0]*beta[0]-36*(lambda*B[1]^2-lambda*B[2]^2-mu^2/lambda)*a[1]*alpha[1]^2-36*a[1]*beta[0]^2 = 0

 (1)
> 

LENGTH := 330:  # to adjust to the page width in your latex file

printlevel:=2:
P := NULL:
L := [op(lhs(eq3))]:
k := 0:
while L <> []  do
  k := k+1:
  l := 0:
  i := 1:
  while length(l) < LENGTH and i <= numelems(L) do
    l := l + L[i]:
    i := i+1:
  end do:
  P := P, l:
  L := L[i..-1]:
end do:

P := [P]:
for i from 1 to numelems(P) do
  p := P[i]:
  if substring(convert(p, string), 1..1) <> "-" then
    P[i] := `#mo("1",mathcolor="white")` + p
  end if:
end do:

P[-1] := P[-1] = 0:
print~(P):

-(9/4)*mu^2*alpha[1]^2*a[1]+3*beta[0]^2*alpha[0]*a[2]+10*beta[0]^2*alpha[0]^3*a[4]+6*beta[0]^2*alpha[0]^2*a[3]-(1/4)*lambda*beta[0]^2*a[1]+(1/4)*(-6*(lambda*B[1]^2-lambda*B[2]^2-mu^2/lambda)*lambda+12*mu^2)*alpha[1]^2*a[1]-w*beta[0]^2

  

`#mo("1",mathcolor="white")`+4*(lambda*B[1]^2-lambda*B[2]^2-mu^2/lambda)*alpha[1]^2*lambda*a[5]*alpha[0]-20*mu*beta[0]*lambda*alpha[1]^4*a[4]+24*mu*beta[0]*alpha[1]^2*alpha[0]*a[3]-30*lambda*beta[0]^2*alpha[1]^2*alpha[0]*a[4]+60*mu*beta[0]*alpha[1]^2*alpha[0]^2*a[4]

  

-7*mu*beta[0]*lambda*a[5]*alpha[1]^2-(lambda*B[1]^2-lambda*B[2]^2-mu^2/lambda)*w*alpha[1]^2+(-2*(lambda*B[1]^2-lambda*B[2]^2-mu^2/lambda)*lambda+4*mu^2)*alpha[1]^4*a[3]+4*(-2*(lambda*B[1]^2-lambda*B[2]^2-mu^2/lambda)*lambda+4*mu^2)*alpha[1]^2*a[5]*alpha[0]

  

-12*mu^2*alpha[1]^2*a[5]*alpha[0]+3*(lambda*B[1]^2-lambda*B[2]^2-mu^2/lambda)*alpha[1]^2*alpha[0]*a[2]+(1/2)*(lambda*B[1]^2-lambda*B[2]^2-mu^2/lambda)*alpha[1]^2*lambda*a[1]+5*(-2*(lambda*B[1]^2-lambda*B[2]^2-mu^2/lambda)*lambda+4*mu^2)*alpha[1]^4*alpha[0]*a[4]

  

`#mo("1",mathcolor="white")`+10*(lambda*B[1]^2-lambda*B[2]^2-mu^2/lambda)*alpha[1]^2*alpha[0]^3*a[4]+6*(lambda*B[1]^2-lambda*B[2]^2-mu^2/lambda)*alpha[1]^2*alpha[0]^2*a[3]-6*lambda*beta[0]^2*alpha[1]^2*a[3]-2*lambda*beta[0]^2*a[5]*alpha[0]+6*mu*beta[0]*alpha[1]^2*a[2]

  

`#mo("1",mathcolor="white")`+3*mu*beta[0]*a[5]*alpha[0]^2+(3/2)*mu*a[1]*alpha[0]*beta[0]-36*(lambda*B[1]^2-lambda*B[2]^2-mu^2/lambda)*a[1]*alpha[1]^2-36*a[1]*beta[0]^2 = 0

 (2)
> 

LENGTH := 150:  # to adjust to the page width in your latex file

printlevel:=2:
P := NULL:
L := [op(lhs(eq3))]:
k := 0:
while L <> [] and k < 20 do
  k := k+1:
  l := 0:
  i := 1:
  while length(l) < LENGTH and i <= numelems(L) do
    l := l + L[i]:
    i := i+1:
  end do:
  P := P, l:
  L := L[i..-1]
end do:

P := [P]:
for i from 1 to numelems(P) do
  p := P[i]:
  if substring(convert(p, string), 1..1) <> "-" then
    P[i] := `#mo("1",mathcolor="white")` + p
  end if:
end do:

P[-1] := P[-1] = 0:
print~(P):

-(9/4)*mu^2*alpha[1]^2*a[1]+3*beta[0]^2*alpha[0]*a[2]+10*beta[0]^2*alpha[0]^3*a[4]+6*beta[0]^2*alpha[0]^2*a[3]

  

-(1/4)*lambda*beta[0]^2*a[1]+(1/4)*(-6*(lambda*B[1]^2-lambda*B[2]^2-mu^2/lambda)*lambda+12*mu^2)*alpha[1]^2*a[1]

  

-w*beta[0]^2+4*(lambda*B[1]^2-lambda*B[2]^2-mu^2/lambda)*alpha[1]^2*lambda*a[5]*alpha[0]-20*mu*beta[0]*lambda*alpha[1]^4*a[4]

  

-30*lambda*a[4]*alpha[0]*alpha[1]^2*beta[0]^2+60*mu*a[4]*alpha[0]^2*alpha[1]^2*beta[0]+24*mu*a[3]*alpha[0]*alpha[1]^2*beta[0]

  

-7*mu*beta[0]*lambda*a[5]*alpha[1]^2-(lambda*B[1]^2-lambda*B[2]^2-mu^2/lambda)*w*alpha[1]^2+(-2*(lambda*B[1]^2-lambda*B[2]^2-mu^2/lambda)*lambda+4*mu^2)*alpha[1]^4*a[3]

  

`#mo("1",mathcolor="white")`+4*(-2*(lambda*B[1]^2-lambda*B[2]^2-mu^2/lambda)*lambda+4*mu^2)*alpha[1]^2*a[5]*alpha[0]-12*mu^2*alpha[1]^2*a[5]*alpha[0]

  

`#mo("1",mathcolor="white")`+3*(lambda*B[1]^2-lambda*B[2]^2-mu^2/lambda)*alpha[1]^2*alpha[0]*a[2]+(1/2)*(lambda*B[1]^2-lambda*B[2]^2-mu^2/lambda)*alpha[1]^2*lambda*a[1]

  

`#mo("1",mathcolor="white")`+5*(-2*(lambda*B[1]^2-lambda*B[2]^2-mu^2/lambda)*lambda+4*mu^2)*alpha[1]^4*alpha[0]*a[4]+10*(lambda*B[1]^2-lambda*B[2]^2-mu^2/lambda)*alpha[1]^2*alpha[0]^3*a[4]

  

`#mo("1",mathcolor="white")`+6*(lambda*B[1]^2-lambda*B[2]^2-mu^2/lambda)*alpha[1]^2*alpha[0]^2*a[3]-6*lambda*beta[0]^2*alpha[1]^2*a[3]

  

-2*lambda*beta[0]^2*a[5]*alpha[0]+6*mu*beta[0]*alpha[1]^2*a[2]+3*mu*beta[0]*a[5]*alpha[0]^2+(3/2)*mu*a[1]*alpha[0]*beta[0]

  

-36*(lambda*B[1]^2-lambda*B[2]^2-mu^2/lambda)*a[1]*alpha[1]^2-36*a[1]*beta[0]^2 = 0

 (3)
> 

# LaTeX output

cat(map(p -> cat(latex(p, output=string), "\\\\"), P)[]);


# If you have to remove the \mbox {{\tt `\#mo(\"1\",mathcolor=\"white\")`}} text:


StringTools:-SubstituteAll(%, "\\mbox {{\\tt `\\#mo(""1"",mathcolor=""white"")`}}", ""):

"-9/4\,{\mu}^{2}{\alpha_{{1}}}^{2}a_{{1}}+3\,{\beta_{{0}}}^{2}\alpha_{{0}}a_{{2}}+10\,{\beta_{{0}}}^{2}{\alpha_{{0}}}^{3}a_{{4}}+6\,{\beta_{{0}}}^{2}{\alpha_{{0}}}^{2}a_{{3}}\\-1/4\,\lambda\,{\beta_{{0}}}^{2}a_{{1}}+1/4\, \left( -6\, \left( \lambda\,{B_{{1}}}^{2}-\lambda\,{B_{{2}}}^{2}-{\frac {{\mu}^{2}}{\lambda}} \right) \lambda+12\,{\mu}^{2} \right) {\alpha_{{1}}}^{2}a_{{1}}\\-w{\beta_{{0}}}^{2}+4\, \left( \lambda\,{B_{{1}}}^{2}-\lambda\,{B_{{2}}}^{2}-{\frac {{\mu}^{2}}{\lambda}} \right) {\alpha_{{1}}}^{2}\lambda\,a_{{5}}\alpha_{{0}}-20\,\mu\,\beta_{{0}}\lambda\,{\alpha_{{1}}}^{4}a_{{4}}\\-30\,\lambda\,{\beta_{{0}}}^{2}{\alpha_{{1}}}^{2}\alpha_{{0}}a_{{4}}+60\,\mu\,\beta_{{0}}{\alpha_{{1}}}^{2}{\alpha_{{0}}}^{2}a_{{4}}+24\,\mu\,\beta_{{0}}{\alpha_{{1}}}^{2}\alpha_{{0}}a_{{3}}\\-7\,\mu\,\beta_{{0}}\lambda\,a_{{5}}{\alpha_{{1}}}^{2}- \left( \lambda\,{B_{{1}}}^{2}-\lambda\,{B_{{2}}}^{2}-{\frac {{\mu}^{2}}{\lambda}} \right) w{\alpha_{{1}}}^{2}+ \left( -2\, \left( \lambda\,{B_{{1}}}^{2}-\lambda\,{B_{{2}}}^{2}-{\frac {{\mu}^{2}}{\lambda}} \right) \lambda+4\,{\mu}^{2} \right) {\alpha_{{1}}}^{4}a_{{3}}\\\mbox {{\tt `\#mo("1",mathcolor="white")`}}+4\, \left( -2\, \left( \lambda\,{B_{{1}}}^{2}-\lambda\,{B_{{2}}}^{2}-{\frac {{\mu}^{2}}{\lambda}} \right) \lambda+4\,{\mu}^{2} \right) {\alpha_{{1}}}^{2}a_{{5}}\alpha_{{0}}-12\,{\mu}^{2}{\alpha_{{1}}}^{2}a_{{5}}\alpha_{{0}}\\\mbox {{\tt `\#mo("1",mathcolor="white")`}}+3\, \left( \lambda\,{B_{{1}}}^{2}-\lambda\,{B_{{2}}}^{2}-{\frac {{\mu}^{2}}{\lambda}} \right) {\alpha_{{1}}}^{2}\alpha_{{0}}a_{{2}}+1/2\, \left( \lambda\,{B_{{1}}}^{2}-\lambda\,{B_{{2}}}^{2}-{\frac {{\mu}^{2}}{\lambda}} \right) {\alpha_{{1}}}^{2}\lambda\,a_{{1}}\\\mbox {{\tt `\#mo("1",mathcolor="white")`}}+5\, \left( -2\, \left( \lambda\,{B_{{1}}}^{2}-\lambda\,{B_{{2}}}^{2}-{\frac {{\mu}^{2}}{\lambda}} \right) \lambda+4\,{\mu}^{2} \right) {\alpha_{{1}}}^{4}\alpha_{{0}}a_{{4}}+10\, \left( \lambda\,{B_{{1}}}^{2}-\lambda\,{B_{{2}}}^{2}-{\frac {{\mu}^{2}}{\lambda}} \right) {\alpha_{{1}}}^{2}{\alpha_{{0}}}^{3}a_{{4}}\\\mbox {{\tt `\#mo("1",mathcolor="white")`}}+6\, \left( \lambda\,{B_{{1}}}^{2}-\lambda\,{B_{{2}}}^{2}-{\frac {{\mu}^{2}}{\lambda}} \right) {\alpha_{{1}}}^{2}{\alpha_{{0}}}^{2}a_{{3}}-6\,\lambda\,{\beta_{{0}}}^{2}{\alpha_{{1}}}^{2}a_{{3}}\\-2\,\lambda\,{\beta_{{0}}}^{2}a_{{5}}\alpha_{{0}}+6\,\mu\,\beta_{{0}}{\alpha_{{1}}}^{2}a_{{2}}+3\,\mu\,\beta_{{0}}a_{{5}}{\alpha_{{0}}}^{2}+3/2\,\mu\,a_{{1}}\alpha_{{0}}\beta_{{0}}\\-36\, \left( \lambda\,{B_{{1}}}^{2}-\lambda\,{B_{{2}}}^{2}-{\frac {{\mu}^{2}}{\lambda}} \right) a_{{1}}{\alpha_{{1}}}^{2}-36\,a_{{1}}{\beta_{{0}}}^{2}=0\\"

 (4)
 

 

Download split.mw

A procedure to emulate operator D...

mmcdara 7284
December 15 2024
1 15

(and nothing more)

restart
alias(F=F(x, t), G=G(x, t)):
                             
with(PDEtools):
undeclare(prime):

ND := proc(F, G, U) 
  local v, w, f, g, a:
  v := op(F):
  if v[1] in U then w := -v[1] else w := v[1] end if:
  if v[2] in U then w := w, -v[2] else w := w, v[2] end if:
  f := op(0, F):
  g := op(0, G):
  a := diff(f(w)*g(v), U);
  convert(subs([w]=~[v], a), diff)
end proc:



An_idea.mw

By the way: I just sent you HERE a reply to your Better simplification and latex export question (LaTeX issue only)

1 2 3 4 5 6 7 Last Page 3 of 61