Advanced Search

dharr

Dr. David Harrington

9339 Reputation

22 Badges

21 years, 335 days
University of Victoria
Professor or university staff
Victoria, British Columbia, Canada

Social Networks and Content at Maplesoft.com

Maple Application Center
Contact Dr. David Harrington
I am a retired professor of chemistry at the University of Victoria, BC, Canada. My research areas are electrochemistry and surface science. I have been a user of Maple since about 1990.

MaplePrimes Activity


These are answers submitted by dharr

new variable names...

dharr 9339
Yesterday at 2:49 PM
1 1

If I understand you correctly, you want to temporarily replace the names, use Ranking or other commands and then convert back to the indexed form. Under the assumption you know what the indexed names are, the following is a simple solution. It has the disadvantage that cat produces global names, so must avoid any global names you are already using.

> 

restart;
> 

odes := {seq(diff(a[i](t),t)=-i*a[i](t), i = 1..3)};

{diff(a[1](t), t) = -a[1](t), diff(a[2](t), t) = -2*a[2](t), diff(a[3](t), t) = -3*a[3](t)}

> 

odes2 := subs(seq(a[i]=cat(qqa__,i),i=1..3),odes);

{diff(qqa__1(t), t) = -qqa__1(t), diff(qqa__2(t), t) = -2*qqa__2(t), diff(qqa__3(t), t) = -3*qqa__3(t)}

> 

odes3:=subs(seq(cat(qqa__,i)=a[i],i=1..3),odes2);

{diff(a[1](t), t) = -a[1](t), diff(a[2](t), t) = -2*a[2](t), diff(a[3](t), t) = -3*a[3](t)}

> 

NULL

Download cat.mw

ways...

dharr 9339
Yesterday at 1:53 PM
3 0
> 

restart;
> 

eqs:={f(1) = x,
      f(2) = 2*x + 2*sin(x),        
      f(n) = 2*f(n-1) - f(n-2) + 2*sin((n-1)*x)/(n-1)};

{f(1) = x, f(2) = 2*x+2*sin(x), f(n) = 2*f(n-1)-f(n-2)+2*sin((n-1)*x)/(n-1)}

rsolve can give a general formula

> 

rsolve(eqs, f(n)) assuming n::posint;

-4*sin(x)-x+(2*sin(x)+x)*(n+1)+csc((1/2)*x)*((-(1/2)*n*LerchPhi(exp(-I*x), 1, n)+1/2)*exp(-((1/2)*I)*x*(-1+2*n))+(-(1/2)*n*LerchPhi(exp(I*x), 1, n)+1/2)*exp(((1/2)*I)*x*(-1+2*n))+(1/2)*exp(-((1/2)*I)*x*(1+2*n))*n*LerchPhi(exp(-I*x), 1, n)+(1/2)*exp(((1/2)*I)*x*(1+2*n))*n*LerchPhi(exp(I*x), 1, n)-I*ln(1-exp(-I*x))*sin((1/2)*x)*n+I*ln(-exp(I*x)+1)*sin((1/2)*x)*n+(n-1)*cos((3/2)*x)-n*cos((1/2)*x))

But you can also ask rsolve to make a procedure

> 

F := rsolve(eqs, f(n), 'makeproc'):

So we can use seq

> 

seq(F(i), i=3..7)

3*x+4*sin(x)+sin(2*x), 4*x+6*sin(x)+2*sin(2*x)+(2/3)*sin(3*x), 5*x+8*sin(x)+3*sin(2*x)+(4/3)*sin(3*x)+(1/2)*sin(4*x), 6*x+10*sin(x)+4*sin(2*x)+2*sin(3*x)+sin(4*x)+(2/5)*sin(5*x), 7*x+12*sin(x)+5*sin(2*x)+(8/3)*sin(3*x)+(3/2)*sin(4*x)+(4/5)*sin(5*x)+(1/3)*sin(6*x)

Or we can use a for loop

> 

for i from 3 to 7 do
  i,F(i);
end do;

3, 3*x+4*sin(x)+sin(2*x)

4, 4*x+6*sin(x)+2*sin(2*x)+(2/3)*sin(3*x)

5, 5*x+8*sin(x)+3*sin(2*x)+(4/3)*sin(3*x)+(1/2)*sin(4*x)

6, 6*x+10*sin(x)+4*sin(2*x)+2*sin(3*x)+sin(4*x)+(2/5)*sin(5*x)

7, 7*x+12*sin(x)+5*sin(2*x)+(8/3)*sin(3*x)+(3/2)*sin(4*x)+(4/5)*sin(5*x)+(1/3)*sin(6*x)

Or we can program this procedure recursively

> 

G := proc(n) option remember;
  if n = 1 then x
  elif n= 2 then 2*x + 2*sin(x)
  else 2*thisproc(n-1) - thisproc(n-2) + 2*sin((n-1)*x)/(n-1)
  end if
end proc:

Try it out

> 

for i from 3 to 7 do
  G(i);
end do;

3*x+4*sin(x)+sin(2*x)

4*x+6*sin(x)+2*sin(2*x)+(2/3)*sin(3*x)

5*x+8*sin(x)+3*sin(2*x)+(4/3)*sin(3*x)+(1/2)*sin(4*x)

6*x+10*sin(x)+4*sin(2*x)+2*sin(3*x)+sin(4*x)+(2/5)*sin(5*x)

7*x+12*sin(x)+5*sin(2*x)+(8/3)*sin(3*x)+(3/2)*sin(4*x)+(4/5)*sin(5*x)+(1/3)*sin(6*x)

> 

NULL

Download recurse.mw

differential quadrature...

dharr 9339
July 23 2026
3 3

My recollection is that Maplesim uses the Method of Lines, suitable for elliptic, parabolic and hyperbolic pdes to convert the pde to a system of odes. This leverages the power of the ode solvers, and is fairly general in terms of the systems it solves. This is probably the best known method of this type.

A similar idea is used in the relatively obscure differential quadrature method, which also produces a system of odes, and can handle nonlinear pdes. For some reason I no longer remember (ease of coding perhaps?), I chose this method and implemented it in Maple a long time ago, and give the code below.

It has some disadvantages: (1) it is clunky in how you enter the pde, (2) the matrix system it solves internally is Vandermonde and so has numerical issues, which could be improved by use of orthogonal polynomials, (3) the output in not directly in a nice form for plotting, (4) there is no error assessment or control - increase the number of space points to check convergence. I see I didn't do much argument checking. The method can be extended to pde systems. Anyway, for what it is worth, here it is:

> 

restart

Procedure to numerically solve a partial differential equation with initial and boundary conditions by the differential quadrature method.

     de     -  a partial differential equation (expression = 0 or equation), with first variable as a "time" variable (in the ode system to be generated), and second variable the space variable. This is written in a special Matrix/Vector form,

                 where D2 is a Matrix operator meaning differentiate wrt the space variable, Dt (treat as Vector) means the 1st derivative of the time variable, Dtt (Vector) means the second derivative,

                  and &* means pointwise multiplication of vectors

     ini      -  list of initial condition(s) - first is a function of the spatial variable that gives the value of dependent variable at time zero, second (if present) is the derivatives of the dependent varable with respect to time at time zero.

     bcs   -   list of the boundary conditions at first and last x values. Here use [1] on the dep. variable to indicate its value at the first x value, and [-1] to indicate the last value.

     vars  -  name of the dependent variable as a function of the independent variables (time and then space)

     grid   - values of the space variable to use - a Vector of values to evaluate the solution at, with the first and last being the boundaries. Need not be equally spaced.

     params - set of parameter values to use

 

Exmple from G. Naadimuthu et al, J. Math. Anal. App., 98 (1984), 220-235. doi: 10.1016/0022-247X(84)90290-7:

 

c[t] = c[x, x]/Npc-c[x]-r*c^2 with boundary conditions `and`(0 = Pe-c(t, 0)+c[x](t, 0)/Npc, 0 = c[x](t, 48)) and initial condition c(0, x) = 0 solved for some x values from 0 to 48. Call as

 

dqsolve(Dt-D2/Npc.D2.c+D2.c+`&*`(r*c, c), [proc (x) options operator, arrow; 0 end proc], [Pe-c[1]+D2/Npc.c, D2.c], c(t, x), Vector([0., 2., 5., 6., 10., 20., 30., 40., 48.]), {Npc = 2, Pe = 0.7e-1, r = 1})

                

> 

dqsolve := proc (de, ini, bcs, vars, grid, params) local C, B, A, U, Ut, Utt, N, t, x, u, des, ins, eq1, eqN, i; N := LinearAlgebra:-Dimension(grid); C := Matrix(N, N, proc (i, j) options operator, arrow; grid[j]^(i-1) end proc); B := Matrix(N, N, proc (i, j) options operator, arrow; if i = 1 then 0 else (i-1)*grid[j]^(i-2) end if end proc); A := LinearAlgebra:-Transpose(LinearAlgebra:-LinearSolve(C, B)); u, t, x := op(0 .. 2, vars); U := Vector(N, proc (i) options operator, arrow; u[i](t) end proc); Ut := Vector(N, proc (i) options operator, arrow; diff(u[i](t), t) end proc); Utt := Vector(N, proc (i) options operator, arrow; diff(u[i](t), t, t) end proc); eval(de, {D2 = A, Dt = Ut, Dtt = Utt, x = grid, u = U}); eval(%, `&*` = `~`[`*`]); des := convert(%, list); ins := seq(u[i](0) = ini[1](grid[i]), i = 1 .. N); if 1 < nops(ini) then ins := ins, seq((D(u[i]))(0) = ini[2](grid[i]), i = 2 .. N-1) end if; eval(bcs[1], {D2 = LinearAlgebra:-Row(A, 1), Dt = Ut, Dtt = Utt, x = grid, u = U}); eq1 := eval(%, `&*` = `~`[`*`]); eval(bcs[2], {D2 = LinearAlgebra:-Row(A, N), Dt = Ut, Dtt = Utt, x = grid, u = U}); eqN := eval(%, `&*` = `~`[`*`]); dsolve(subs(params, {op(des[2 .. N-1]), ins, eq1, eqN}), convert(U, set), numeric, method = rosenbrock_dae) end proc
> 

ans1 := dqsolve(Dt-D2/Npc.D2.c+D2.c+`&*`(r*c, c), [proc (x) options operator, arrow; 0 end proc], [Pe-c[1]+D2/Npc.c, D2.c], c(t, x), Vector([0., 2., 5., 6., 10., 20., 30., 40., 48.]), {Npc = 2, Pe = 0.7e-1, r = 1})

proc (x_rosenbrock_dae) local _res, _dat, _vars, _solnproc, _xout, _ndsol, _pars, _n, _i; option `Copyright - do not display`; if 1 < nargs then error "invalid input: too many arguments" end if; _EnvDSNumericSaveDigits := Digits; Digits := 15; if _EnvInFsolve = true then _xout := evalf[_EnvDSNumericSaveDigits](x_rosenbrock_dae) else _xout := evalf(x_rosenbrock_dae) end if; _dat := Array(1..4, {(1) = proc (_x_in) local _x_out, _dtbl, _dat, _vmap, _x0, _y0, _val, _digits, _neq, _nevar, _ndisc, _nevt, _pars, _ini, _par, _i, _j, _k, _src, _t1, _tmap; option `Copyright (c) 2002 by Waterloo Maple Inc. All rights reserved.`; table( [( "complex" ) = false ] ) _x_out := _x_in; _pars := []; _dtbl := array( 1 .. 4, [( 1 ) = (array( 1 .. 29, [( 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..65, {(1) = 9, (2) = 7, (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) = 2, (23) = 3, (24) = 0, (25) = 1, (26) = 15, (27) = 1, (28) = 0, (29) = 1, (30) = 3, (31) = 3, (32) = 0, (33) = 2, (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, (55) = 0, (56) = 0, (57) = 0, (58) = 0, (59) = 10000, (60) = 0, (61) = 1000, (62) = 0, (63) = 0, (64) = -1, (65) = 0}, datatype = integer[8])), ( 5 ) = (Array(1..28, {(1) = .0, (2) = 0.10e-5, (3) = .0, (4) = 0.500001e-14, (5) = .0, (6) = 0.6974672852195059e-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..9, {(1) = .0, (2) = .0, (3) = .0, (4) = .0, (5) = .0, (6) = .0, (7) = .0, (8) = 0.45241164318878924e-1, (9) = -2.50772573876595}, 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..9, {(1) = 1.0, (2) = 1.0, (3) = 1.0, (4) = 1.0, (5) = 1.0, (6) = 1.0, (7) = 1.0, (8) = 1.0, (9) = 1.0}, datatype = float[8], order = C_order), Array(1..7, {(1) = .0, (2) = .0, (3) = .0, (4) = .0, (5) = .0, (6) = .0, (7) = .0}, datatype = float[8], order = C_order), Array(1..7, {(1) = .0, (2) = .0, (3) = .0, (4) = .0, (5) = .0, (6) = .0, (7) = .0}, datatype = float[8], order = C_order), Array(1..7, {(1) = .0, (2) = .0, (3) = .0, (4) = .0, (5) = .0, (6) = .0, (7) = .0}, datatype = float[8], order = C_order), Array(1..7, {(1) = .0, (2) = .0, (3) = .0, (4) = .0, (5) = .0, (6) = .0, (7) = .0}, datatype = float[8], order = C_order), Array(1..7, 1..7, {(1, 1) = .3329759696413988, (1, 2) = -1.580399771016325, (1, 3) = 1.1815730217272935, (1, 4) = -0.8344732629988748e-1, (1, 5) = 0.33150451583014275e-2, (1, 6) = -0.4751169720728332e-3, (1, 7) = 0.5661622332835021e-4, (2, 1) = .16429516335725994, (2, 2) = .5895708905984295, (2, 3) = -.7755144052164096, (2, 4) = 0.3883199070687317e-1, (2, 5) = -0.11899774556985534e-2, (2, 6) = 0.15940983553975295e-3, (2, 7) = -0.184406988633596e-4, (3, 1) = -0.15770649761358037e-1, (3, 2) = 1.0858447827454358, (3, 3) = -1.0790050234168924, (3, 4) = 0.7785786648967076e-2, (3, 5) = -0.5046092870994577e-4, (3, 6) = 0.35223424731640414e-5, (3, 7) = -0.26812476503604694e-6, (4, 1) = .15737122084608388, (4, 2) = -1.238811317007148, (4, 3) = 1.5339825115208023, (4, 4) = -.4273054035850302, (4, 5) = -0.37068790352788004e-2, (4, 6) = 0.45470806968293093e-3, (4, 7) = -0.50562734998277276e-4, (5, 1) = -1.1081167700501322, (5, 2) = 5.963143313030375, (5, 3) = -5.959208312587414, (5, 4) = 1.1026453023888734, (5, 5) = -.16729830416704958, (5, 6) = -0.1629559718584923e-1, (5, 7) = 0.14153129416528889e-2, (6, 1) = 3.163469061842351, (6, 2) = -15.784596246809281, (6, 3) = 15.286761497856096, (6, 4) = -2.373938897126878, (6, 5) = .42753328829811055, (6, 6) = -.1573908620582749, (6, 7) = -0.12320597040825701e-1, (7, 1) = -12.33708238191245, (7, 2) = 59.60067182886593, (7, 3) = -57.00625011660849, (7, 4) = 8.330759017294627, (7, 5) = -1.0911724204165487, (7, 6) = .5156639724151103, (7, 7) = -.20281854337134586}, datatype = float[8], order = C_order), Array(1..7, 1..7, {(1, 1) = .0, (1, 2) = .0, (1, 3) = .0, (1, 4) = .0, (1, 5) = .0, (1, 6) = .0, (1, 7) = .0, (2, 1) = .0, (2, 2) = .0, (2, 3) = .0, (2, 4) = .0, (2, 5) = .0, (2, 6) = .0, (2, 7) = .0, (3, 1) = .0, (3, 2) = .0, (3, 3) = .0, (3, 4) = .0, (3, 5) = .0, (3, 6) = .0, (3, 7) = .0, (4, 1) = .0, (4, 2) = .0, (4, 3) = .0, (4, 4) = .0, (4, 5) = .0, (4, 6) = .0, (4, 7) = .0, (5, 1) = .0, (5, 2) = .0, (5, 3) = .0, (5, 4) = .0, (5, 5) = .0, (5, 6) = .0, (5, 7) = .0, (6, 1) = .0, (6, 2) = .0, (6, 3) = .0, (6, 4) = .0, (6, 5) = .0, (6, 6) = .0, (6, 7) = .0, (7, 1) = .0, (7, 2) = .0, (7, 3) = .0, (7, 4) = .0, (7, 5) = .0, (7, 6) = .0, (7, 7) = .0}, datatype = float[8], order = C_order), Array(1..7, {(1) = .0, (2) = .0, (3) = .0, (4) = .0, (5) = .0, (6) = .0, (7) = .0}, datatype = float[8], order = C_order), Array(1..7, 1..7, {(1, 1) = .0, (1, 2) = .0, (1, 3) = .0, (1, 4) = .0, (1, 5) = .0, (1, 6) = .0, (1, 7) = .0, (2, 1) = .0, (2, 2) = .0, (2, 3) = .0, (2, 4) = .0, (2, 5) = .0, (2, 6) = .0, (2, 7) = .0, (3, 1) = .0, (3, 2) = .0, (3, 3) = .0, (3, 4) = .0, (3, 5) = .0, (3, 6) = .0, (3, 7) = .0, (4, 1) = .0, (4, 2) = .0, (4, 3) = .0, (4, 4) = .0, (4, 5) = .0, (4, 6) = .0, (4, 7) = .0, (5, 1) = .0, (5, 2) = .0, (5, 3) = .0, (5, 4) = .0, (5, 5) = .0, (5, 6) = .0, (5, 7) = .0, (6, 1) = .0, (6, 2) = .0, (6, 3) = .0, (6, 4) = .0, (6, 5) = .0, (6, 6) = .0, (6, 7) = .0, (7, 1) = .0, (7, 2) = .0, (7, 3) = .0, (7, 4) = .0, (7, 5) = .0, (7, 6) = .0, (7, 7) = .0}, datatype = float[8], order = C_order), Array(1..9, 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, (7, 1) = .0, (7, 2) = .0, (7, 3) = .0, (7, 4) = .0, (7, 5) = .0, (7, 6) = .0, (8, 1) = .0, (8, 2) = .0, (8, 3) = .0, (8, 4) = .0, (8, 5) = .0, (8, 6) = .0, (9, 1) = .0, (9, 2) = .0, (9, 3) = .0, (9, 4) = .0, (9, 5) = .0, (9, 6) = .0}, datatype = float[8], order = C_order), Array(1..7, {(1) = 0, (2) = 0, (3) = 0, (4) = 0, (5) = 0, (6) = 0, (7) = 0}, datatype = integer[8]), Array(1..9, {(1) = .0, (2) = .0, (3) = .0, (4) = .0, (5) = .0, (6) = .0, (7) = .0, (8) = .0, (9) = .0}, datatype = float[8], order = C_order), Array(1..9, {(1) = .0, (2) = .0, (3) = .0, (4) = .0, (5) = .0, (6) = .0, (7) = .0, (8) = .0, (9) = .0}, datatype = float[8], order = C_order), Array(1..9, {(1) = .0, (2) = .0, (3) = .0, (4) = .0, (5) = .0, (6) = .0, (7) = .0, (8) = .0, (9) = .0}, datatype = float[8], order = C_order), Array(1..9, {(1) = .0, (2) = .0, (3) = .0, (4) = .0, (5) = .0, (6) = .0, (7) = .0, (8) = .0, (9) = .0}, datatype = float[8], order = C_order), Array(1..7, {(1) = .0, (2) = .0, (3) = .0, (4) = .0, (5) = .0, (6) = .0, (7) = .0}, datatype = float[8], order = C_order), Array(1..9, {(1) = .0, (2) = .0, (3) = .0, (4) = .0, (5) = .0, (6) = .0, (7) = .0, (8) = .0, (9) = .0}, datatype = float[8], order = C_order), Array(1..18, {(1) = .0, (2) = .0, (3) = .0, (4) = .0, (5) = .0, (6) = .0, (7) = .0, (8) = .0, (9) = .0, (10) = .0, (11) = .0, (12) = .0, (13) = .0, (14) = .0, (15) = .0, (16) = .0, (17) = .0, (18) = .0}, datatype = float[8], order = C_order), Array(1..7, {(1) = 0, (2) = 0, (3) = 0, (4) = 0, (5) = 0, (6) = 0, (7) = 0}, datatype = integer[8])]), ( 8 ) = ([Array(1..9, {(1) = .0, (2) = .0, (3) = .0, (4) = .0, (5) = .0, (6) = .0, (7) = .0, (8) = undefined, (9) = undefined}, datatype = float[8], order = C_order), Array(1..9, {(1) = .0, (2) = .0, (3) = .0, (4) = .0, (5) = .0, (6) = .0, (7) = .0, (8) = .0, (9) = .0}, datatype = float[8], order = C_order), Array(1..7, {(1) = 0.5493644742035918e-2, (2) = 0.8809292579333025e-4, (3) = -0.9381855610225828e-4, (4) = 0.11946658948398413e-2, (5) = -0.4993557334618059e-2, (6) = 0.13748590981476145e-1, (7) = -0.6578127618609739e-1}, datatype = float[8], order = C_order), 0, 0]), ( 11 ) = (Array(1..6, 0..9, {(1, 1) = .0, (1, 2) = .0, (1, 3) = .0, (1, 4) = .0, (1, 5) = .0, (1, 6) = .0, (1, 7) = .0, (1, 8) = .0, (1, 9) = .0, (2, 0) = .0, (2, 1) = .0, (2, 2) = .0, (2, 3) = .0, (2, 4) = .0, (2, 5) = .0, (2, 6) = .0, (2, 7) = .0, (2, 8) = .0, (2, 9) = .0, (3, 0) = .0, (3, 1) = .0, (3, 2) = .0, (3, 3) = .0, (3, 4) = .0, (3, 5) = .0, (3, 6) = .0, (3, 7) = .0, (3, 8) = .0, (3, 9) = .0, (4, 0) = .0, (4, 1) = .0, (4, 2) = .0, (4, 3) = .0, (4, 4) = .0, (4, 5) = .0, (4, 6) = .0, (4, 7) = .0, (4, 8) = .0, (4, 9) = .0, (5, 0) = .0, (5, 1) = .0, (5, 2) = .0, (5, 3) = .0, (5, 4) = .0, (5, 5) = .0, (5, 6) = .0, (5, 7) = .0, (5, 8) = .0, (5, 9) = .0, (6, 0) = .0, (6, 1) = .0, (6, 2) = .0, (6, 3) = .0, (6, 4) = .0, (6, 5) = .0, (6, 6) = .0, (6, 7) = .0, (6, 8) = .0, (6, 9) = .0}, datatype = float[8], order = C_order)), ( 10 ) = ([proc (N, X, Y, YP) Y[8] := (3514832466093085443090617069166623/5338067473797507403186260903279717)*Y[1]-(5592630595056987680362164098594656/5338067473797507403186260903279717)*Y[2]+(4275583989866991789668571457329134/5338067473797507403186260903279717)*Y[3]-(65938250012025170785072859619179/1094988199753334851935643262211224)*Y[4]+(286952739535255523150903656086907/113878772774346824601306899269967296)*Y[5]-(251077430790343497377557306236571/683272636646080947607841395619803776)*Y[6]+(40221952585783356328985763030053/911030182194774596810455194159738368)*Y[7]+6037509693183398696508499523796992/133451686844937685079656522581992925; Y[9] := (355612602727572159400502484847165440/1779355824599169134395420301093239)*Y[1]-(1694529111474084569310722945670578176/1779355824599169134395420301093239)*Y[2]+(1612477683508971505328721123784392704/1779355824599169134395420301093239)*Y[3]-(17701300773866934552102962032721920/136873524969166856491955407776403)*Y[4]+(27356353014537693701834188451857920/1779355824599169134395420301093239)*Y[5]-(9983303898426828356330533555796992/1779355824599169134395420301093239)*Y[6]+(4706847717305859813379651926367432/1779355824599169134395420301093239)*Y[7]-111553409994261185467478379623088128/44483895614979228359885507527330975; if N < 1 then return 0 end if; YP[1] := -Y[1]^2+(10673220026814313437841540788775729405289328700171/32054024914497359015055286024514817969982774706176)*Y[1]-(12664543408755801016665888159661234293688408224845/8013506228624339753763821506128704492495693676544)*Y[2]+(2367135692296537174953610542368314597178284658673/2003376557156084938440955376532176123123923419136)*Y[3]-(823022361927625883004839722700127688964938956617/9862776896768418158478549546004559375379315294208)*Y[4]+(6800674566200179935542566347253080997316704544413/2051457594527830976963538305568948350078897581195264)*Y[5]-(15594917130366097178024526053694981335646338169713/32823321512445295631416612889103173601262361299124224)*Y[6]+(3716665002253689430759187407859306759547586226229/65646643024890591262833225778206347202524722598248448)*Y[7]+455875824692730076882593070682643/44483895614979228359885507527330975; YP[2] := -Y[2]^2+(14173590012078339510635964418156178811004576921445/24040518685873019261291464518386113477487081029632)*Y[2]-(9684559803299787878414898970342925138575113927375/525173144199124730102665806225650777620197780785987584)*Y[7]+(31394161705004039218011623467365846148003552150349/196939929074671773788499677334619041607574167794745344)*Y[6]-(39059012620955250372934594865631941804218776766339/32823321512445295631416612889103173601262361299124224)*Y[5]+(12255720345576788653914862085043365308936465890367/315608860696589381071313585472145900012138089414656)*Y[4]+(252783420460081103674923479345313056828237192687261/1538593195895873232722653729176711262559173185896448)*Y[1]-(596600593592601390656552644012024108084731380862443/769296597947936616361326864588355631279586592948224)*Y[3]-602894247350213222091482772292819/533806747379750740318626090327971700; YP[3] := -Y[3]^2+(104417087158197518775268884981464934329804076433669/96162074743492077045165858073544453909948324118528)*Y[2]-(35202981472908075718308352322847523214865785555/131293286049781182525666451556412694405049445196496896)*Y[7]+(1387379753683260424705747362031034113471082077871/393879858149343547576999354669238083215148335589490688)*Y[6]-(207036910857891427909338070356419858499097471655/4102915189055661953927076611137896700157795162390528)*Y[5]+(1228631626953607562203217365307732002137211875805/157804430348294690535656792736072950006069044707328)*Y[4]-(3033076802210294116474227400074449428639500795127/192324149486984154090331716147088907819896648237056)*Y[1]-(415037446841674473228573000012727915363329415175333/384648298973968308180663432294177815639793296474112)*Y[3]+178209748394230546028843271334549/2135226989519002961274504361311886800; YP[4] := -Y[4]^2-(9927222204927101390975859501148936001166866716767/8013506228624339753763821506128704492495693676544)*Y[2]-(3319273814794050928095263221276793649190408472167/65646643024890591262833225778206347202524722598248448)*Y[7]+(14925029165506223697202905269380683960054124325135/32823321512445295631416612889103173601262361299124224)*Y[6]-(60836041191349557837348557990209044926960790329019/16411660756222647815708306444551586800631180649562112)*Y[5]-(67430685797483848620760530830612363868539237274789/157804430348294690535656792736072950006069044707328)*Y[4]+(40355048270601908598955137440870431313605360917615/256432199315978872120442288196118543759862197649408)*Y[1]+(196681254570764133889111602733779230653068695327523/128216099657989436060221144098059271879931098824704)*Y[3]-273202198146864530819410581441387/177935582459916913439542030109323900; YP[5] := -Y[5]^2+(11946421520287122380162659053785240965719104639013/2003376557156084938440955376532176123123923419136)*Y[2]+(725863620696835959309093418639495763955717319369/512864398631957744240884576392237087519724395298816)*Y[7]-(33429726604276753545913990204371322460555039067851/2051457594527830976963538305568948350078897581195264)*Y[6]-(171602688317560448016961229690326989933050710597575/1025728797263915488481769152784474175039448790597632)*Y[5]+(10875144613731206503200078061770936018267852633247/9862776896768418158478549546004559375379315294208)*Y[4]-(17759801277679640902517745223059063742894720366521/16027012457248679507527643012257408984991387353088)*Y[1]-(47754152930589191754448076464211054909510285773553/8013506228624339753763821506128704492495693676544)*Y[3]+572065295026737768152651088217604/44483895614979228359885507527330975; YP[6] := -Y[6]^2-(7905622521257910519988003129519124321905125107413/500844139289021234610238844133044030780980854784)*Y[2]-(25275182368519006693166587386159367979776174615215/2051457594527830976963538305568948350078897581195264)*Y[7]-(80720169819682584888868208145745997700810103102943/512864398631957744240884576392237087519724395298816)*Y[6]+(109633301399076950016728181869249393236099884113077/256432199315978872120442288196118543759862197649408)*Y[5]-(5853407427230717142701734183694420865810506895661/2465694224192104539619637386501139843844828823552)*Y[4]+(12675239515567039944551297519881934381406370083851/4006753114312169876881910753064352246247846838272)*Y[1]+(30625139619641147960405273274042028359510343573647/2003376557156084938440955376532176123123923419136)*Y[3]-1711126743444247862808574537168032/44483895614979228359885507527330975; YP[7] := -Y[7]^2+(44775970774763648179348695950577401852286957412247/751266208933531851915358266199566046171471282176)*Y[2]-(104018410277554926782470739232745389535109542107835/512864398631957744240884576392237087519724395298816)*Y[7]+(198349269831631489529667977174558397252241764416573/384648298973968308180663432294177815639793296474112)*Y[6]-(69952935900088878647833542450704994948969563804445/64108049828994718030110572049029635939965549412352)*Y[5]+(5135276098019913448088068173376248456898525211917/616423556048026134904909346625284960961207205888)*Y[4]-(9268433110360033470880850866535782836890040616819/751266208933531851915358266199566046171471282176)*Y[1]-(85653738821242330865355571813265752983593027690897/1502532417867063703830716532399132092342942564352)*Y[3]+20460279495760496627856440737517824/133451686844937685079656522581992925; 0 end proc, proc (X, Y, FX, FY) FX[1 .. 7] := 0; FY[1 .. 7, 1 .. 7] := 0; FY[1, 1] := -2*Y[1]+10673220026814313437841540788775729405289328700171/32054024914497359015055286024514817969982774706176; FY[2, 1] := 252783420460081103674923479345313056828237192687261/1538593195895873232722653729176711262559173185896448; FY[3, 1] := -3033076802210294116474227400074449428639500795127/192324149486984154090331716147088907819896648237056; FY[4, 1] := 40355048270601908598955137440870431313605360917615/256432199315978872120442288196118543759862197649408; FY[5, 1] := -17759801277679640902517745223059063742894720366521/16027012457248679507527643012257408984991387353088; FY[6, 1] := 12675239515567039944551297519881934381406370083851/4006753114312169876881910753064352246247846838272; FY[7, 1] := -9268433110360033470880850866535782836890040616819/751266208933531851915358266199566046171471282176; FY[1, 2] := -12664543408755801016665888159661234293688408224845/8013506228624339753763821506128704492495693676544; FY[2, 2] := -2*Y[2]+14173590012078339510635964418156178811004576921445/24040518685873019261291464518386113477487081029632; FY[3, 2] := 104417087158197518775268884981464934329804076433669/96162074743492077045165858073544453909948324118528; FY[4, 2] := -9927222204927101390975859501148936001166866716767/8013506228624339753763821506128704492495693676544; FY[5, 2] := 11946421520287122380162659053785240965719104639013/2003376557156084938440955376532176123123923419136; FY[6, 2] := -7905622521257910519988003129519124321905125107413/500844139289021234610238844133044030780980854784; FY[7, 2] := 44775970774763648179348695950577401852286957412247/751266208933531851915358266199566046171471282176; FY[1, 3] := 2367135692296537174953610542368314597178284658673/2003376557156084938440955376532176123123923419136; FY[2, 3] := -596600593592601390656552644012024108084731380862443/769296597947936616361326864588355631279586592948224; FY[3, 3] := -2*Y[3]-415037446841674473228573000012727915363329415175333/384648298973968308180663432294177815639793296474112; FY[4, 3] := 196681254570764133889111602733779230653068695327523/128216099657989436060221144098059271879931098824704; FY[5, 3] := -47754152930589191754448076464211054909510285773553/8013506228624339753763821506128704492495693676544; FY[6, 3] := 30625139619641147960405273274042028359510343573647/2003376557156084938440955376532176123123923419136; FY[7, 3] := -85653738821242330865355571813265752983593027690897/1502532417867063703830716532399132092342942564352; FY[1, 4] := -823022361927625883004839722700127688964938956617/9862776896768418158478549546004559375379315294208; FY[2, 4] := 12255720345576788653914862085043365308936465890367/315608860696589381071313585472145900012138089414656; FY[3, 4] := 1228631626953607562203217365307732002137211875805/157804430348294690535656792736072950006069044707328; FY[4, 4] := -2*Y[4]-67430685797483848620760530830612363868539237274789/157804430348294690535656792736072950006069044707328; FY[5, 4] := 10875144613731206503200078061770936018267852633247/9862776896768418158478549546004559375379315294208; FY[6, 4] := -5853407427230717142701734183694420865810506895661/2465694224192104539619637386501139843844828823552; FY[7, 4] := 5135276098019913448088068173376248456898525211917/616423556048026134904909346625284960961207205888; FY[1, 5] := 6800674566200179935542566347253080997316704544413/2051457594527830976963538305568948350078897581195264; FY[2, 5] := -39059012620955250372934594865631941804218776766339/32823321512445295631416612889103173601262361299124224; FY[3, 5] := -207036910857891427909338070356419858499097471655/4102915189055661953927076611137896700157795162390528; FY[4, 5] := -60836041191349557837348557990209044926960790329019/16411660756222647815708306444551586800631180649562112; FY[5, 5] := -2*Y[5]-171602688317560448016961229690326989933050710597575/1025728797263915488481769152784474175039448790597632; FY[6, 5] := 109633301399076950016728181869249393236099884113077/256432199315978872120442288196118543759862197649408; FY[7, 5] := -69952935900088878647833542450704994948969563804445/64108049828994718030110572049029635939965549412352; FY[1, 6] := -15594917130366097178024526053694981335646338169713/32823321512445295631416612889103173601262361299124224; FY[2, 6] := 31394161705004039218011623467365846148003552150349/196939929074671773788499677334619041607574167794745344; FY[3, 6] := 1387379753683260424705747362031034113471082077871/393879858149343547576999354669238083215148335589490688; FY[4, 6] := 14925029165506223697202905269380683960054124325135/32823321512445295631416612889103173601262361299124224; FY[5, 6] := -33429726604276753545913990204371322460555039067851/2051457594527830976963538305568948350078897581195264; FY[6, 6] := -2*Y[6]-80720169819682584888868208145745997700810103102943/512864398631957744240884576392237087519724395298816; FY[7, 6] := 198349269831631489529667977174558397252241764416573/384648298973968308180663432294177815639793296474112; FY[1, 7] := 3716665002253689430759187407859306759547586226229/65646643024890591262833225778206347202524722598248448; FY[2, 7] := -9684559803299787878414898970342925138575113927375/525173144199124730102665806225650777620197780785987584; FY[3, 7] := -35202981472908075718308352322847523214865785555/131293286049781182525666451556412694405049445196496896; FY[4, 7] := -3319273814794050928095263221276793649190408472167/65646643024890591262833225778206347202524722598248448; FY[5, 7] := 725863620696835959309093418639495763955717319369/512864398631957744240884576392237087519724395298816; FY[6, 7] := -25275182368519006693166587386159367979776174615215/2051457594527830976963538305568948350078897581195264; FY[7, 7] := -2*Y[7]-104018410277554926782470739232745389535109542107835/512864398631957744240884576392237087519724395298816; 0 end proc, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]), ( 13 ) = (), ( 12 ) = (), ( 15 ) = ("rosenbrock"), ( 14 ) = ([0, 0]), ( 18 ) = ([]), ( 19 ) = (0), ( 16 ) = ([0, 0, 0, 0, 0, 0, []]), ( 17 ) = ([proc (N, X, Y, YP) Y[8] := (3514832466093085443090617069166623/5338067473797507403186260903279717)*Y[1]-(5592630595056987680362164098594656/5338067473797507403186260903279717)*Y[2]+(4275583989866991789668571457329134/5338067473797507403186260903279717)*Y[3]-(65938250012025170785072859619179/1094988199753334851935643262211224)*Y[4]+(286952739535255523150903656086907/113878772774346824601306899269967296)*Y[5]-(251077430790343497377557306236571/683272636646080947607841395619803776)*Y[6]+(40221952585783356328985763030053/911030182194774596810455194159738368)*Y[7]+6037509693183398696508499523796992/133451686844937685079656522581992925; Y[9] := (355612602727572159400502484847165440/1779355824599169134395420301093239)*Y[1]-(1694529111474084569310722945670578176/1779355824599169134395420301093239)*Y[2]+(1612477683508971505328721123784392704/1779355824599169134395420301093239)*Y[3]-(17701300773866934552102962032721920/136873524969166856491955407776403)*Y[4]+(27356353014537693701834188451857920/1779355824599169134395420301093239)*Y[5]-(9983303898426828356330533555796992/1779355824599169134395420301093239)*Y[6]+(4706847717305859813379651926367432/1779355824599169134395420301093239)*Y[7]-111553409994261185467478379623088128/44483895614979228359885507527330975; if N < 1 then return 0 end if; YP[1] := -Y[1]^2+(10673220026814313437841540788775729405289328700171/32054024914497359015055286024514817969982774706176)*Y[1]-(12664543408755801016665888159661234293688408224845/8013506228624339753763821506128704492495693676544)*Y[2]+(2367135692296537174953610542368314597178284658673/2003376557156084938440955376532176123123923419136)*Y[3]-(823022361927625883004839722700127688964938956617/9862776896768418158478549546004559375379315294208)*Y[4]+(6800674566200179935542566347253080997316704544413/2051457594527830976963538305568948350078897581195264)*Y[5]-(15594917130366097178024526053694981335646338169713/32823321512445295631416612889103173601262361299124224)*Y[6]+(3716665002253689430759187407859306759547586226229/65646643024890591262833225778206347202524722598248448)*Y[7]+455875824692730076882593070682643/44483895614979228359885507527330975; YP[2] := -Y[2]^2+(14173590012078339510635964418156178811004576921445/24040518685873019261291464518386113477487081029632)*Y[2]-(9684559803299787878414898970342925138575113927375/525173144199124730102665806225650777620197780785987584)*Y[7]+(31394161705004039218011623467365846148003552150349/196939929074671773788499677334619041607574167794745344)*Y[6]-(39059012620955250372934594865631941804218776766339/32823321512445295631416612889103173601262361299124224)*Y[5]+(12255720345576788653914862085043365308936465890367/315608860696589381071313585472145900012138089414656)*Y[4]+(252783420460081103674923479345313056828237192687261/1538593195895873232722653729176711262559173185896448)*Y[1]-(596600593592601390656552644012024108084731380862443/769296597947936616361326864588355631279586592948224)*Y[3]-602894247350213222091482772292819/533806747379750740318626090327971700; YP[3] := -Y[3]^2+(104417087158197518775268884981464934329804076433669/96162074743492077045165858073544453909948324118528)*Y[2]-(35202981472908075718308352322847523214865785555/131293286049781182525666451556412694405049445196496896)*Y[7]+(1387379753683260424705747362031034113471082077871/393879858149343547576999354669238083215148335589490688)*Y[6]-(207036910857891427909338070356419858499097471655/4102915189055661953927076611137896700157795162390528)*Y[5]+(1228631626953607562203217365307732002137211875805/157804430348294690535656792736072950006069044707328)*Y[4]-(3033076802210294116474227400074449428639500795127/192324149486984154090331716147088907819896648237056)*Y[1]-(415037446841674473228573000012727915363329415175333/384648298973968308180663432294177815639793296474112)*Y[3]+178209748394230546028843271334549/2135226989519002961274504361311886800; YP[4] := -Y[4]^2-(9927222204927101390975859501148936001166866716767/8013506228624339753763821506128704492495693676544)*Y[2]-(3319273814794050928095263221276793649190408472167/65646643024890591262833225778206347202524722598248448)*Y[7]+(14925029165506223697202905269380683960054124325135/32823321512445295631416612889103173601262361299124224)*Y[6]-(60836041191349557837348557990209044926960790329019/16411660756222647815708306444551586800631180649562112)*Y[5]-(67430685797483848620760530830612363868539237274789/157804430348294690535656792736072950006069044707328)*Y[4]+(40355048270601908598955137440870431313605360917615/256432199315978872120442288196118543759862197649408)*Y[1]+(196681254570764133889111602733779230653068695327523/128216099657989436060221144098059271879931098824704)*Y[3]-273202198146864530819410581441387/177935582459916913439542030109323900; YP[5] := -Y[5]^2+(11946421520287122380162659053785240965719104639013/2003376557156084938440955376532176123123923419136)*Y[2]+(725863620696835959309093418639495763955717319369/512864398631957744240884576392237087519724395298816)*Y[7]-(33429726604276753545913990204371322460555039067851/2051457594527830976963538305568948350078897581195264)*Y[6]-(171602688317560448016961229690326989933050710597575/1025728797263915488481769152784474175039448790597632)*Y[5]+(10875144613731206503200078061770936018267852633247/9862776896768418158478549546004559375379315294208)*Y[4]-(17759801277679640902517745223059063742894720366521/16027012457248679507527643012257408984991387353088)*Y[1]-(47754152930589191754448076464211054909510285773553/8013506228624339753763821506128704492495693676544)*Y[3]+572065295026737768152651088217604/44483895614979228359885507527330975; YP[6] := -Y[6]^2-(7905622521257910519988003129519124321905125107413/500844139289021234610238844133044030780980854784)*Y[2]-(25275182368519006693166587386159367979776174615215/2051457594527830976963538305568948350078897581195264)*Y[7]-(80720169819682584888868208145745997700810103102943/512864398631957744240884576392237087519724395298816)*Y[6]+(109633301399076950016728181869249393236099884113077/256432199315978872120442288196118543759862197649408)*Y[5]-(5853407427230717142701734183694420865810506895661/2465694224192104539619637386501139843844828823552)*Y[4]+(12675239515567039944551297519881934381406370083851/4006753114312169876881910753064352246247846838272)*Y[1]+(30625139619641147960405273274042028359510343573647/2003376557156084938440955376532176123123923419136)*Y[3]-1711126743444247862808574537168032/44483895614979228359885507527330975; YP[7] := -Y[7]^2+(44775970774763648179348695950577401852286957412247/751266208933531851915358266199566046171471282176)*Y[2]-(104018410277554926782470739232745389535109542107835/512864398631957744240884576392237087519724395298816)*Y[7]+(198349269831631489529667977174558397252241764416573/384648298973968308180663432294177815639793296474112)*Y[6]-(69952935900088878647833542450704994948969563804445/64108049828994718030110572049029635939965549412352)*Y[5]+(5135276098019913448088068173376248456898525211917/616423556048026134904909346625284960961207205888)*Y[4]-(9268433110360033470880850866535782836890040616819/751266208933531851915358266199566046171471282176)*Y[1]-(85653738821242330865355571813265752983593027690897/1502532417867063703830716532399132092342942564352)*Y[3]+20460279495760496627856440737517824/133451686844937685079656522581992925; 0 end proc, proc (X, Y, FX, FY) FX[1 .. 7] := 0; FY[1 .. 7, 1 .. 7] := 0; FY[1, 1] := -2*Y[1]+10673220026814313437841540788775729405289328700171/32054024914497359015055286024514817969982774706176; FY[2, 1] := 252783420460081103674923479345313056828237192687261/1538593195895873232722653729176711262559173185896448; FY[3, 1] := -3033076802210294116474227400074449428639500795127/192324149486984154090331716147088907819896648237056; FY[4, 1] := 40355048270601908598955137440870431313605360917615/256432199315978872120442288196118543759862197649408; FY[5, 1] := -17759801277679640902517745223059063742894720366521/16027012457248679507527643012257408984991387353088; FY[6, 1] := 12675239515567039944551297519881934381406370083851/4006753114312169876881910753064352246247846838272; FY[7, 1] := -9268433110360033470880850866535782836890040616819/751266208933531851915358266199566046171471282176; FY[1, 2] := -12664543408755801016665888159661234293688408224845/8013506228624339753763821506128704492495693676544; FY[2, 2] := -2*Y[2]+14173590012078339510635964418156178811004576921445/24040518685873019261291464518386113477487081029632; FY[3, 2] := 104417087158197518775268884981464934329804076433669/96162074743492077045165858073544453909948324118528; FY[4, 2] := -9927222204927101390975859501148936001166866716767/8013506228624339753763821506128704492495693676544; FY[5, 2] := 11946421520287122380162659053785240965719104639013/2003376557156084938440955376532176123123923419136; FY[6, 2] := -7905622521257910519988003129519124321905125107413/500844139289021234610238844133044030780980854784; FY[7, 2] := 44775970774763648179348695950577401852286957412247/751266208933531851915358266199566046171471282176; FY[1, 3] := 2367135692296537174953610542368314597178284658673/2003376557156084938440955376532176123123923419136; FY[2, 3] := -596600593592601390656552644012024108084731380862443/769296597947936616361326864588355631279586592948224; FY[3, 3] := -2*Y[3]-415037446841674473228573000012727915363329415175333/384648298973968308180663432294177815639793296474112; FY[4, 3] := 196681254570764133889111602733779230653068695327523/128216099657989436060221144098059271879931098824704; FY[5, 3] := -47754152930589191754448076464211054909510285773553/8013506228624339753763821506128704492495693676544; FY[6, 3] := 30625139619641147960405273274042028359510343573647/2003376557156084938440955376532176123123923419136; FY[7, 3] := -85653738821242330865355571813265752983593027690897/1502532417867063703830716532399132092342942564352; FY[1, 4] := -823022361927625883004839722700127688964938956617/9862776896768418158478549546004559375379315294208; FY[2, 4] := 12255720345576788653914862085043365308936465890367/315608860696589381071313585472145900012138089414656; FY[3, 4] := 1228631626953607562203217365307732002137211875805/157804430348294690535656792736072950006069044707328; FY[4, 4] := -2*Y[4]-67430685797483848620760530830612363868539237274789/157804430348294690535656792736072950006069044707328; FY[5, 4] := 10875144613731206503200078061770936018267852633247/9862776896768418158478549546004559375379315294208; FY[6, 4] := -5853407427230717142701734183694420865810506895661/2465694224192104539619637386501139843844828823552; FY[7, 4] := 5135276098019913448088068173376248456898525211917/616423556048026134904909346625284960961207205888; FY[1, 5] := 6800674566200179935542566347253080997316704544413/2051457594527830976963538305568948350078897581195264; FY[2, 5] := -39059012620955250372934594865631941804218776766339/32823321512445295631416612889103173601262361299124224; FY[3, 5] := -207036910857891427909338070356419858499097471655/4102915189055661953927076611137896700157795162390528; FY[4, 5] := -60836041191349557837348557990209044926960790329019/16411660756222647815708306444551586800631180649562112; FY[5, 5] := -2*Y[5]-171602688317560448016961229690326989933050710597575/1025728797263915488481769152784474175039448790597632; FY[6, 5] := 109633301399076950016728181869249393236099884113077/256432199315978872120442288196118543759862197649408; FY[7, 5] := -69952935900088878647833542450704994948969563804445/64108049828994718030110572049029635939965549412352; FY[1, 6] := -15594917130366097178024526053694981335646338169713/32823321512445295631416612889103173601262361299124224; FY[2, 6] := 31394161705004039218011623467365846148003552150349/196939929074671773788499677334619041607574167794745344; FY[3, 6] := 1387379753683260424705747362031034113471082077871/393879858149343547576999354669238083215148335589490688; FY[4, 6] := 14925029165506223697202905269380683960054124325135/32823321512445295631416612889103173601262361299124224; FY[5, 6] := -33429726604276753545913990204371322460555039067851/2051457594527830976963538305568948350078897581195264; FY[6, 6] := -2*Y[6]-80720169819682584888868208145745997700810103102943/512864398631957744240884576392237087519724395298816; FY[7, 6] := 198349269831631489529667977174558397252241764416573/384648298973968308180663432294177815639793296474112; FY[1, 7] := 3716665002253689430759187407859306759547586226229/65646643024890591262833225778206347202524722598248448; FY[2, 7] := -9684559803299787878414898970342925138575113927375/525173144199124730102665806225650777620197780785987584; FY[3, 7] := -35202981472908075718308352322847523214865785555/131293286049781182525666451556412694405049445196496896; FY[4, 7] := -3319273814794050928095263221276793649190408472167/65646643024890591262833225778206347202524722598248448; FY[5, 7] := 725863620696835959309093418639495763955717319369/512864398631957744240884576392237087519724395298816; FY[6, 7] := -25275182368519006693166587386159367979776174615215/2051457594527830976963538305568948350078897581195264; FY[7, 7] := -2*Y[7]-104018410277554926782470739232745389535109542107835/512864398631957744240884576392237087519724395298816; 0 end proc, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]), ( 22 ) = (0), ( 23 ) = (0), ( 20 ) = ([]), ( 21 ) = (0), ( 27 ) = (""), ( 26 ) = (Array(1..0, {})), ( 25 ) = (Array(1..0, {})), ( 24 ) = (0), ( 29 ) = (0), ( 28 ) = (0) ] )) ] ); _y0 := Array(0..9, {(1) = 0., (2) = 0., (3) = 0., (4) = 0., (5) = 0., (6) = 0., (7) = 0., (8) = 0., (9) = Float(undefined)}); _vmap := array( 1 .. 9, [( 1 ) = (8), ( 2 ) = (1), ( 3 ) = (2), ( 4 ) = (3), ( 5 ) = (4), ( 6 ) = (5), ( 7 ) = (6), ( 9 ) = (9), ( 8 ) = (7) ] ); _x0 := _dtbl[1][5][5]; _neq := _dtbl[1][4][1]; _nevar := _dtbl[1][4][3]; _ndisc := _dtbl[1][4][4]; _nevt := _dtbl[1][4][16]; _tmap := _dtbl[1][29]; if not type(_x_out, 'numeric') then if member(_x_out, ["start", "left", "right"]) then if _Env_smart_dsolve_numeric = true or _dtbl[1][4][10] = 1 then if _x_out = "left" then if type(_dtbl[2], 'array') then return _dtbl[2][5][1] end if elif _x_out = "right" then if type(_dtbl[3], 'array') then return _dtbl[3][5][1] end if end if end if; return _dtbl[1][5][5] elif _x_out = "method" then return _dtbl[1][15] elif _x_out = "storage" then return evalb(_dtbl[1][4][10] = 1) elif _x_out = "leftdata" then if not type(_dtbl[2], 'array') then return NULL else return eval(_dtbl[2]) end if elif _x_out = "rightdata" then if not type(_dtbl[3], 'array') then return NULL else return eval(_dtbl[3]) end if elif _x_out = "enginedata" then return _dtbl[1] elif _x_out = "enginereset" then _dtbl[2] := evaln(_dtbl[2]); _dtbl[3] := evaln(_dtbl[3]); return NULL elif _x_out = "initial" then return procname(_y0[0]) elif _x_out = "laxtol" then return _dtbl[`if`(member(_dtbl[4], {2, 3}), _dtbl[4], 1)][5][18] elif _x_out = "numfun" then return `if`(member(_dtbl[4], {2, 3}), _dtbl[_dtbl[4]][4][18], 0) elif _x_out = "parameters" then if _tmap = 0 then return [seq(_y0[_neq+_i], _i = 1 .. nops(_pars))] else return [seq(`dsolve/numeric/FloatToDiscrete`(_y0[_neq+_i], _tmap[_neq+_i]), _i = 1 .. nops(_pars))] end if elif _x_out = "initial_and_parameters" then if _tmap = 0 then return procname(_y0[0]), [seq(_y0[_neq+_i], _i = 1 .. nops(_pars))] else return procname(_y0[0]), [seq(`dsolve/numeric/FloatToDiscrete`(_y0[_neq+_i], _tmap[_neq+_i]), _i = 1 .. nops(_pars))] end if elif _x_out = "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 _x_out := _dtbl[_dtbl[4]][5][1] end if elif _x_out = "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 _x_out = "map" then return copy(_vmap) elif type(_x_in, `=`) and type(rhs(_x_in), 'list') and member(lhs(_x_in), {"initial", "parameters", "initial_and_parameters"}) then _ini, _par := [], []; if lhs(_x_in) = "initial" then _ini := rhs(_x_in) elif lhs(_x_in) = "parameters" then _par := rhs(_x_in) elif select(type, rhs(_x_in), `=`) <> [] then _par, _ini := selectremove(type, rhs(_x_in), `=`) elif nops(rhs(_x_in)) < nops(_pars)+1 then error "insufficient data for specification of initial and parameters" else _par := rhs(_x_in)[-nops(_pars) .. -1]; _ini := rhs(_x_in)[1 .. -nops(_pars)-1] end if; _x_out := lhs(_x_out); _i := false; if _par <> [] then _i := `dsolve/numeric/process_parameters`(_neq, _pars, _par, _y0, _tmap) end if; if _ini <> [] then _i := `dsolve/numeric/process_initial`(_neq-_nevar, _ini, _y0, _pars, _vmap) or _i end if; if _i then `dsolve/numeric/SC/reinitialize`(_dtbl, _y0, _neq, 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 end if; if _x_out = "initial" then return [_y0[0], seq(_y0[_vmap[_i]], _i = 1 .. _neq-_nevar)] elif _x_out = "parameters" then return [seq(_y0[_neq+_i], _i = 1 .. nops(_pars))] else return [_y0[0], seq(_y0[_vmap[_i]], _i = 1 .. _neq-_nevar)], [seq(_y0[_neq+_i], _i = 1 .. nops(_pars))] end if elif _x_in = "eventstop" then if _nevt = 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 _x_in = "eventstatus" then if _nevt = 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][_nevt+1, 1]))})]; return ':-enabled' = _i[1], ':-disabled' = _i[2] elif _x_in = "eventclear" then if _nevt = 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 _nevt < _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 _nevt 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][_nevt+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(_x_in, `=`) and member(lhs(_x_in), {"eventdisable", "eventenable"}) then if _nevt = 0 then error "this solution has no events" end if; if type(rhs(_x_in), {('list')('posint'), ('set')('posint')}) then _i := {op(rhs(_x_in))} elif type(rhs(_x_in), 'posint') then _i := {rhs(_x_in)} else error "event identifiers must be integers in the range 1..%1", round(_dtbl[1][3][1][_nevt+1, 1]) end if; if select(proc (a) options operator, arrow; _nevt < a end proc, _i) <> {} then error "event identifiers must be integers in the range 1..%1", round(_dtbl[1][3][1][_nevt+1, 1]) end if; _k := {}; for _j to _nevt do if member(round(_dtbl[1][3][1][_j, 1]), _i) then _k := `union`(_k, {_j}) end if end do; _i := _k; if lhs(_x_in) = "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 _nevt+1 do if _k <= _nevt 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 _nevt+1 do if _k <= _nevt 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(_x_in, `=`) and lhs(_x_in) = "eventfired" then if not type(rhs(_x_in), 'list') then error "'eventfired' must be specified as a list" end if; if _nevt = 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(_x_in) 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][_nevt+1, 1]) < _src then error "event identifiers must be integers in the range 1..%1", round(_dtbl[1][3][1][_nevt+1, 1]) end if; _src := {seq(`if`(_dtbl[1][3][1][_j, 1]-_src = 0., _j, NULL), _j = 1 .. _nevt)}; 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(_x_in, `=`) and lhs(_x_in) = "direction" then if not member(rhs(_x_in), {-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(_x_in), {1, ':-right'}), 3, 2); _dtbl[4] := _i; _dtbl[_i] := `dsolve/numeric/SC/IVPdcopy`(_dtbl[1], `if`(assigned(_dtbl[_i]), _dtbl[_i], NULL)); if 0 < _nevt then for _j to _nevt+1 do if _j <= _nevt 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 _x_in = "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][_nevt+1, 12]) end if elif type(_x_in, `=`) and lhs(_x_in) = "setdatacallback" then if not type(rhs(_x_in), 'nonegint') then error "data callback must be a nonnegative integer (address)" end if; _dtbl[1][28] := rhs(_x_in) elif type(_x_in, `=`) and lhs(_x_in) = "Array" then _i := rhs(_x_in); if not type(_i, 'list') or nops(_i) <> 2 or not type(_i[1], 'numeric') or not type(_i[2], 'posint') or _i[2] < 2 then error "Array output must be specified as [end time, min number of points]" end if; _src := array(1 .. 1, [`dsolve/numeric/SC/IVPdcopy`(_dtbl[1])]); if 0 < 0 then if `dsolve/numeric/checkglobals`(op(_src[1][14]), _pars, _neq, _y0) then `dsolve/numeric/SC/reinitialize`(_src, _y0, _neq, procname, _pars, 1) end if end if; if _src[1][4][7] = 0 then error "parameters must be initialized before solution can be computed" end if; _val := `dsolve/numeric/SC/IVPvalues`(_src[1], _x0 .. _i[1], _i[2], _i[2], []); if _val[3] <> "" then `dsolve/numeric/warning`(cat(`requested integration incomplete, received error:`, convert(_val[3], 'symbol'))) end if; _t1 := Array(1 .. _val[2], 0 .. _neq-_nevar+nops(_pars)); _t1[() .. (), 0] := _val[1][1 .. _val[2], 0]; for _i to _neq-_nevar do _t1[() .. (), _i] := _val[1][1 .. _val[2], _vmap[_i]] end do; for _i to nops(_pars) do _t1[() .. (), _neq-_nevar+_i] := _val[1][1 .. _val[2], _neq+_i] end do; return _t1 else return "procname" end if end if; if _x_out = _x0 then return [_x0, seq(evalf(_dtbl[1][6][_vmap[_i]]), _i = 1 .. _neq-_nevar)] end if; _i := `if`(_x0 <= _x_out, 3, 2); if _x_in = "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 .. _neq-_nevar-_ndisc), seq(_dat[8][1][_vmap[_i]], _i = _neq-_nevar-_ndisc+1 .. _neq-_nevar)] 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 < _nevt then for _j to _nevt+1 do if _j <= _nevt 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 _x_in <> "last" then if 0 < 0 then if `dsolve/numeric/checkglobals`(op(_dtbl[1][14]), _pars, _neq, _y0) then `dsolve/numeric/SC/reinitialize`(_dtbl, _y0, _neq, 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, _x_out) catch: userinfo(2, `dsolve/debug`, print(`Exception in solnproc:`, [lastexception][2 .. -1])); error end try; if _dat[17] <> _dtbl[1][17] then _dtbl[1][17] := _dat[17]; _dtbl[1][10] := _dat[10] end if; if _src = 0 and 100 < _dat[4][9] then _val := _dat[3][1][_nevt+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] := _x_out else _dtbl[1][5][1] := _val end if; if _i = 3 and _val < _x_out 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 = _nevt+1 then error "constraint projection failure on event at t=%1", evalf[8](_val) elif _dat[4][9]-100 = _nevt+2 then error "index-1 and derivative evaluation failure on event at t=%1", evalf[8](_val) elif _dat[4][9]-100 = _nevt+3 then error "maximum number of event iterations reached (%1) at t=%2", round(_dat[3][1][_nevt+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'; _x_out := _val end if else error "cannot evaluate the solution further right of %1", evalf[8](_val) end if elif _i = 2 and _x_out < _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 = _nevt+1 then error "constraint projection failure on event at t=%1", evalf[8](_val) elif _dat[4][9]-100 = _nevt+2 then error "index-1 and derivative evaluation failure on event at t=%1", evalf[8](_val) elif _dat[4][9]-100 = _nevt+3 then error "maximum number of event iterations reached (%1) at t=%2", round(_dat[3][1][_nevt+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'; _x_out := _val end if else error "cannot evaluate the solution further left of %1", evalf[8](_val) end if end if; if _EnvInFsolve = true then _digits := _dat[4][26]; if type(_EnvDSNumericSaveDigits, 'posint') then _dat[4][26] := _EnvDSNumericSaveDigits else _dat[4][26] := Digits end if; _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, _x_out, _src); _dat[4][25] := _i; _dat[4][26] := _digits; [_x_out, seq(_val[_vmap[_i]], _i = 1 .. _neq-_nevar)] else Digits := _dat[4][26]; _val := `dsolve/numeric/SC/IVPval`(eval(_dat, 2), _x_out, _src); [_x_out, seq(_val[_vmap[_i]], _i = 1 .. _neq-_nevar)] end if end proc, (2) = Array(0..0, {}), (3) = [t, c[1](t), c[2](t), c[3](t), c[4](t), c[5](t), c[6](t), c[7](t), c[8](t), c[9](t)], (4) = []}); _vars := _dat[3]; _pars := map(lhs, _dat[4]); _n := nops(_vars)-1; _solnproc := _dat[1]; if not type(_xout, 'numeric') then if member(x_rosenbrock_dae, ["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_rosenbrock_dae, '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_rosenbrock_dae, ["last", 'last', "initial", 'initial', "parameters", 'parameters', "initial_and_parameters", 'initial_and_parameters', NULL]) then _xout := convert(x_rosenbrock_dae, '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_rosenbrock_dae), 'string') = rhs(x_rosenbrock_dae); 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', "Array", 'Array', NULL]) then return _solnproc(convert(lhs(x_rosenbrock_dae), 'string') = rhs(x_rosenbrock_dae)) elif _xout = "solnprocedure" then return eval(_solnproc) elif _xout = "sysvars" then return _vars elif _xout = "Array" then return [op(_vars), op(_pars)] end if; if procname <> unknown then return ('procname')(x_rosenbrock_dae) else _ndsol := 1; _ndsol := _ndsol; _ndsol := pointto(_dat[2][0]); return ('_ndsol')(x_rosenbrock_dae) end if end if; try _res := _solnproc(_xout); [seq(_vars[_i+1] = _res[_i+1], _i = 0 .. _n)] catch: error end try end proc

The values at the chosen grid points in x at time t = 40

> 

ans1(40)

[t = 40., c[1](t) = HFloat(0.06783752347417663), c[2](t) = HFloat(0.0601372829457519), c[3](t) = HFloat(0.051330327808960824), c[4](t) = HFloat(0.04893033513351124), c[5](t) = HFloat(0.0411926398168522), c[6](t) = HFloat(0.029467840801792863), c[7](t) = HFloat(0.022498908869120333), c[8](t) = HFloat(0.014063604681351398), c[9](t) = HFloat(0.0058276425872099856)]

> 

NULL

Download differential_quadrature_Naadimuthu_procedure2.mw

Using Maple with some help...

dharr 9339
July 15 2026
3 0

Here's a step by step way of deriving Mariusz's result (and evaluating the sum).

We want to solveP(n, m) = p*P(n-1, m)+(1-p)*P(n, m-1), P(n, 0) = 0, P(0, m) = 1

> 

restart
> 

rec1 := P(n, m)-p*P(n-1, m)-(1-p)*P(n, m-1)

P(n, m)-p*P(n-1, m)-(1-p)*P(n, m-1)

Make a generating function G(n) = sum(P(n, m)*z^m, m = 0 .. infinity) where z marks m. Do this by multiplying by z^m and summing.

For
m = (0*x^0)*P(n, 0)-p*x^0*P(n-1, 0)-(x*(1-p)*(1/x))*P(n, -1)*m and (0*x^0)*P(n, 0)-p*x^0*P(n-1, 0)-(x*(1-p)*(1/x))*P(n, -1)*m = x*P(n, 1)-p*x*P(n-1, 1)-x*(1-p)*x^0*P(n, 0)

etc,

and summing and assuming P(n, -1) = 0 gives

> 

rec := G(n)-p*G(n-1)-z*(1-p)*G(n)

G(n)-p*G(n-1)-z*(1-p)*G(n)

Now "G(0)=(&sum;)P(0,m)z^(m)=0+z+z^(2)+... =1/(1-z)-1" if we want P(0, 0) = 0. This is the boundary condition for rec, which can be solved by rsolve

> 

gf := simplify(rsolve({rec, G(0) = 1/(1-z)-1}, G(n)))

-z*(p/(1+(p-1)*z))^n/(-1+z)

It will be simpler to solve if we remove the product z/(1-z) which causes partial summation and an offset by 1. (Maple can solve it directly, but it is a mess)

> 

gf2 := simplify((1-z)*gf/z)

(p/(p*z-z+1))^n

To find the coefficient of x^m from the generating function, we differentiate m times, evaluate at z = 0 and divide by factorial(m)

> 

`assuming`([(eval(diff(gf2, `$`(z, m)), z = 0))/factorial(m)], [n::nonnegint, m::nonnegint, p >= 0, p <= 1]); b := convert(%, binomial)

p^n*(p-1)^m*pochhammer(-n-m+1, m)/factorial(m)

p^n*(p-1)^m*binomial(-n, -n-m)

For some reason, Maple likes the negative values in the binomial and (p-1) instead of (1-p). If we look at, say m=5, we see [(n+4)(n+3)...n]/5! or (n+m-1)!/(n-1)!/m! = binomial(n+m-1,m)

> 

(diff(gf2, `$`(z, 5)))/factorial(5); factor(eval(%, z = 0))

-(1/120)*(p/(p*z-z+1))^n*n^5*(p-1)^5/(p*z-z+1)^5-(1/12)*(p/(p*z-z+1))^n*n^4*(p-1)^5/(p*z-z+1)^5-(7/24)*(p/(p*z-z+1))^n*n^3*(p-1)^5/(p*z-z+1)^5-(5/12)*(p/(p*z-z+1))^n*n^2*(p-1)^5/(p*z-z+1)^5-(1/5)*(p/(p*z-z+1))^n*n*(p-1)^5/(p*z-z+1)^5

-(1/120)*p^n*n*(p-1)^5*(n+4)*(n+3)*(n+2)*(n+1)

Suggesting the simpler form:

> 

b2 := p^n*(1-p)^m*binomial(n+m-1, m)

p^n*(1-p)^m*binomial(n+m-1, m)

Though simplify doesn't think so unless we assume the top parts of the binomials are nonnegative (which isn't true for the conditions we want with Maple's binomial)

> 

`assuming`([simplify(b-b2)], [n::nonposint, m::nonnegint, n+m-1 > 0]); simplify([seq(seq(eval(b-b2, {m = i, n = j}), i = 0 .. 4), j = 0 .. 4)])

0

[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]

We removed the partial summation so we need to put it back (summing to m-1 because of the shift).
We can leave it as a summation or we can evaluate the sum in terms of a hypergeometric function.

> 

k := 'k'; P1 := unapply(('add')(eval(b2, m = k), k = 0 .. m-1), n, m); P2 := unapply(sum(eval(b2, m = k), k = 0 .. m-1), n, m)

k

proc (n, m) options operator, arrow; add(p^n*(1-p)^k*binomial(n+k-1, k), k = 0 .. m-1) end proc

proc (n, m) options operator, arrow; 1-p^n*(1-p)^m*binomial(n+m-1, m)*hypergeom([1, n+m], [m+1], 1-p) end proc

Let's check it out

> 

PP := proc(n, m) option remember;
if m<0 or n<0 then 0
elif m=0 then 0 # put m first so P(0,0)=0
elif n=0 then 1
else p*thisproc(n - 1, m) + (1 - p)*thisproc(n, m - 1)
end if;
end proc:

rows indexed by n, columns indexed by m

> 

k := 5; B := expand(Array(0 .. k, 0 .. k, PP)); Tabulate(DataFrame(Matrix(B), columns = [seq(0 .. k)], rows = [seq(0 .. k)]))
> 

expand(PP(4, 5)); expand(P1(4, 5)); expand(simplify(P2(4, 5)))

35*p^8-160*p^7+280*p^6-224*p^5+70*p^4

35*p^8-160*p^7+280*p^6-224*p^5+70*p^4

35*p^8-160*p^7+280*p^6-224*p^5+70*p^4

> 

P2(0, 0)

0

Download rsolve.mw

The DataFrame Table that didn't display:

worksheetdir...

dharr 9339
July 12 2026
1 0

Like this? Of course it can be combined as Import(cat(interface(worksheetdir),"/data.xlsx"));

> 

restart
> 

with(ExcelTools)
> 

dir := interface(worksheetdir)

"C:\Users\dharr\OneDrive - University of Victoria\Desktop"

> 

Import(cat(dir, "/data.xlsx"))

Matrix(%id = 36893490243433306492)

 

NULL

Download testexcel.mw

invisible multiplication...

dharr 9339
July 11 2026
1 2

Here's a hack, by substituting into the MathML

c:=Unit('mm');
subs("&lowast;"=" ",InertForm:-Display(`%*`(Units:-Split(c))));

gives

 

elaboration...

dharr 9339
July 07 2026
5 1

@Rouben Rostamian  Since I just finished a more extended explanation, I'll add it in case it is useful to others

> 

restart;
> 

with(PDETools):

Differential equation over 0 < x and x < 5: The function u(x) is unknown

> 

de := diff(u(x),x) + x*u(x) = 0;

diff(u(x), x)+x*u(x) = 0

Transform the domain from 0 < x and x < 5 to 0 < xi and xi < 1:

> 

tr := x = 5*xi;

x = 5*xi

The result contains u(xi), not u(5*xi) because we don't know the function, it is just some function to be determined.
But implicitly we mean the image of x is the same as the image of the xi value that is equivalent to x.

> 

dchange({tr}, de, {xi});

(1/5)*(diff(u(xi), xi))+5*xi*u(xi) = 0

This is made explicit in the transformation you used. We have a new name for the transformed function, but the result is essentially the same.

> 

tr2 := x = 5*xi, u(x) = v(xi);
dchange({tr2}, de, {xi, v(xi)});

x = 5*xi, u(x) = v(xi)

(1/5)*(diff(v(xi), xi))+5*xi*v(xi) = 0

Now consider a slightly different de, involving the known function sin

> 

de2 := diff(u(x),x) + x*sin(x) = 0;

diff(u(x), x)+x*sin(x) = 0

The transformed equation involves sin(5*xi)

> 

dchange({tr}, de2, {xi});

(1/5)*(diff(u(xi), xi))+5*xi*sin(5*xi) = 0

In the next case, we know the function, we just don't want to tell Maple.

> 

de3 := diff(u(x),x) + x*f(x) = 0;

diff(u(x), x)+x*f(x) = 0

Then the following is not what we want

> 

dchange({tr}, de3, {xi});

(1/5)*(diff(u(xi), xi))+5*xi*f(xi) = 0

But if we tell Maple the function is known we get the required result

> 

dchange({tr}, de3, {xi}, known = f);

(1/5)*(diff(u(xi), xi))+5*xi*f(5*xi) = 0

So a generic first order ODE, if we mean that the function is known

> 

DE := F(x,u(x),diff(u(x),x));

F(x, u(x), diff(u(x), x))

and we tell Maple that, then we get the expected 5*xi

> 

dchange({tr2}, DE, {xi, v(xi)}, known=F);

F(5*xi, v(xi), (1/5)*(diff(v(xi), xi)))

> 

 

Download dchange-problem.mw

local/global confusion...

dharr 9339
July 07 2026
1 1

One solution is to move the ode inside MoveRod as you had in the Rod_triangle code. But a better solution is just to move the ics and dsolve outside MoveRod and then get rid of MoveRod altogether - it was defined and then immediately invoked only once.  Like this:

Rod_parabola.mw

The problem is that when ode is defined outside MoveRod, x is a global. Then inside MoveRod you declare x as a local, so the dsolve command has two x's that look the same but are not. You can confirm this by changing dsolve to ddsolve, then after MoveRod(); add indets(%,function). Towards the end you will see x(t) appears twice.

updates...

dharr 9339
July 06 2026
1 0

You should be able to update to 2025.2 through Tools -> check for updates.

If you mean the Customer Support Updates and the SupportTools package that makes small updates between numbered updates, that was brought in in Maple 2026. I recall it was experimental during 2025, someone may be able to help you install it manually (search Mapleprimes, perhaps?), but probably there wasn't much beyond the 2025.2 update.

series...

dharr 9339
July 06 2026
1 0

The message means there are two variables (nu and x) and it does not know which to expand in. It does not seem to know the general nu case. But this can be done in other ways:

Download Bessel.mw

seq...

dharr 9339
July 04 2026
2 0

Simplest way is seq:

seq(ithprime(b), b = 1..19);

?dsolve, numeric, interactive...

dharr 9339
July 04 2026
2 1

It is hard to find but is at ?dsolve,numeric,interactive.

2016 -> 2026...

dharr 9339
July 01 2026
1 1

I'm not sure I completely understand what you want to do. The simplest is to just open your 2016 worksheet in Maple 2026. You may get a warning message but I think almost everything will work the same; if not you can ask a specific question about that here.

Since you want to use 1D input in 2026, go to the File tab in the top ribbon, and then in the bottom-right corner of the drop-menu choose "options". That brings up a popup menu. Choose its Display tab, and then choose Input Display as Maple Notation and choose Apply Globally or Apply to Session. The output display could also be changed, though that seems more unusual.

If you just want to copy a few lines you can directly copy from 2016 to 2026. If you haven't set your Input Display to Maple Notation, you need to use Ctrl-m on the input line just before you paste to ensure that input line is in 1D mode.

pde2...

dharr 9339
June 21 2026
1 8

As you used, the pde2 paper must be using the ordinary derivative. With some additional assumptions it can give Eq 3.5.

Download pde2.mw

some success...

dharr 9339
June 21 2026
1 0

I'm assuming the notation (lambda - nu)r means pochhammer. If so, straightforward evaluation of the double sum works in some cases.

> 

restart
> 

F := proc (lambda, nu, kappa, mu, alpha, beta, x) options operator, arrow; GAMMA(kappa+lambda+1)*GAMMA(mu+nu+1)*(Sum(Sum(pochhammer(lambda-kappa, r)*pochhammer(nu-mu, s)*BesselJ(kappa+lambda+r, alpha*x)*BesselJ(mu+nu+s, beta*sqrt(1-x^2))*((1/2)*alpha*x)^(kappa+lambda+r)*((1/2)*beta*sqrt(1-x^2))^(mu+nu+s)/(2^(r+s)*factorial(r)*factorial(s)*GAMMA(lambda+r+1)*GAMMA(nu+s+1)), r = 0 .. infinity), s = 0 .. infinity)) end proc

proc (lambda, nu, kappa, mu, alpha, beta, x) options operator, arrow; GAMMA(kappa+lambda+1)*GAMMA(mu+nu+1)*(Sum(Sum(pochhammer(lambda-kappa, r)*pochhammer(nu-mu, s)*BesselJ(kappa+lambda+r, alpha*x)*BesselJ(mu+nu+s, beta*sqrt(1-x^2))*((1/2)*alpha*x)^(kappa+lambda+r)*((1/2)*beta*sqrt(1-x^2))^(mu+nu+s)/(2^(r+s)*factorial(r)*factorial(s)*GAMMA(lambda+r+1)*GAMMA(nu+s+1)), r = 0 .. infinity), s = 0 .. infinity)) end proc

Success:

> 

F(1/2, 1/2, 1/2, 1/2, alpha, beta, x); value(%)

Sum(Sum(pochhammer(0, r)*pochhammer(0, s)*BesselJ(1+r, alpha*x)*BesselJ(1+s, beta*(-x^2+1)^(1/2))*((1/2)*alpha*x)^(1+r)*((1/2)*beta*(-x^2+1)^(1/2))^(1+s)/(2^(r+s)*factorial(r)*factorial(s)*GAMMA(3/2+r)*GAMMA(3/2+s)), r = 0 .. infinity), s = 0 .. infinity)

BesselJ(1, alpha*x)*BesselJ(1, beta*(-x^2+1)^(1/2))*alpha*x*beta*(-x^2+1)^(1/2)/Pi

Success

> 

F(-1/2, 1/2, 1/2, 1/2, alpha, beta, x); value(%)

Sum(Sum(pochhammer(-1, r)*pochhammer(0, s)*BesselJ(r, alpha*x)*BesselJ(1+s, beta*(-x^2+1)^(1/2))*((1/2)*alpha*x)^r*((1/2)*beta*(-x^2+1)^(1/2))^(1+s)/(2^(r+s)*factorial(r)*factorial(s)*GAMMA(1/2+r)*GAMMA(3/2+s)), r = 0 .. infinity), s = 0 .. infinity)

BesselJ(0, alpha*x)*BesselJ(1, beta*(-x^2+1)^(1/2))*beta*(-x^2+1)^(1/2)/Pi-(1/2)*BesselJ(1, alpha*x)*BesselJ(1, beta*(-x^2+1)^(1/2))*alpha*x*beta*(-x^2+1)^(1/2)/Pi

Success

> 

F(1/2, -1/2, 1/2, 1/2, alpha, beta, x); value(%)

Sum(Sum(pochhammer(0, r)*pochhammer(-1, s)*BesselJ(1+r, alpha*x)*BesselJ(s, beta*(-x^2+1)^(1/2))*((1/2)*alpha*x)^(1+r)*((1/2)*beta*(-x^2+1)^(1/2))^s/(2^(r+s)*factorial(r)*factorial(s)*GAMMA(3/2+r)*GAMMA(1/2+s)), r = 0 .. infinity), s = 0 .. infinity)

BesselJ(1, alpha*x)*BesselJ(0, beta*(-x^2+1)^(1/2))*alpha*x/Pi-(1/2)*BesselJ(1, alpha*x)*BesselJ(1, beta*(-x^2+1)^(1/2))*alpha*x*beta*(-x^2+1)^(1/2)/Pi

Reduces to a single sum

> 

F(1/2, 1/2, -1/2, 1/2, alpha, beta, x); value(%)

Sum(Sum(pochhammer(1, r)*pochhammer(0, s)*BesselJ(r, alpha*x)*BesselJ(1+s, beta*(-x^2+1)^(1/2))*((1/2)*alpha*x)^r*((1/2)*beta*(-x^2+1)^(1/2))^(1+s)/(2^(r+s)*factorial(r)*factorial(s)*GAMMA(3/2+r)*GAMMA(3/2+s)), r = 0 .. infinity), s = 0 .. infinity)

sum(pochhammer(1, r)*BesselJ(r, alpha*x)*BesselJ(1, beta*(-x^2+1)^(1/2))*((1/2)*alpha*x)^r*beta*(-x^2+1)^(1/2)/(2^r*factorial(r)*GAMMA(3/2+r)*Pi^(1/2)), r = 0 .. infinity)

> 

NULL

Download Hughes2.mw

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