tomleslie

12482 Reputation

19 Badges

12 years, 226 days

MaplePrimes Activity


These are answers submitted by tomleslie

Maple is not "parsing" your 2-D input in the way that you think. Actual parsing is more akin to performing a symbolic-order differentiation, and symbolic-order differentiation only works for integer values - see more comments in the attached.

The ShowSolution() command is only intended to work for "simple" calculus problems. As its help page states (emphasis added)

The ShowSolution command is used to show the solution steps for a Calculus1 problem, that is, a limit, differentiation or integration problem such as can be expected to be encountered in a single-variable calculus course..

So think of the very first "single-variable calculus course" you ever did, and the nature of the problems which you solved. ShowSolution() will (probably) work for most of them - but you didn't do fractional differentiation in such a course!

f := diff(x, [`$`(x, 1/2)])

x

(1)

#
# OP entered above - presumably using palettes.
# So what has (s)he *actually* entered?
#
  f := diff(x, [x $ 1/2]);
#
# Hmmm so what does the second argument in the
# above evaluate to?
#
  [x $ 1/2];
#
# So OP is actually asking for
#
  diff(x, []);
#
# which Maple interprets as don't differentiate!
#
# This is consistent with the "approved" method
# for performing "symbolic order" differentiation
# which would be
#
  dso:=diff(x, x$n);
#
# which we can now evaluate for any INTEGER n -
# see the help at ?diff, symbolicorder
#  
  eval(dso, n=-3);
  eval(dso, n=-2);
  eval(dso, n=-1);
  eval(dso, n=0);
  eval(dso, n=1);

x

 

[]

 

x

 

pochhammer(2-n, n)*x^(1-n)

 

(1/24)*x^4

 

(1/6)*x^3

 

(1/2)*x^2

 

x

 

1

(2)

f := diff(x, [`$`(x, n)]); eval(dso, n = -3); eval(dso, n = -2); eval(dso, n = -1); eval(dso, n = 0); eval(dso, n = 1)

pochhammer(2-n, n)*x^(1-n)

 

(1/24)*x^4

 

(1/6)*x^3

 

(1/2)*x^2

 

x

 

1

(3)

#
# fracdiff works as it should
#
  fracdiff(x, x, 1/2);

2*x^(1/2)/Pi^(1/2)

(4)

  ShowSolution(fracdiff(x, x, 1/2));
#
# The ShowSolution() command is in the
# Student[Calculus1] and cannot be used in the above
# "short form" unless the Student[Calculus1] package
# has been loaded using a with() command.
#
# However even in "long form", as below
#
  Student[Calculus1]:-ShowSolution(fracdiff(x, x, 1/2));
#
# it still won't work, because it is only intended to
# function with *trivial* calculus problems, of the
# type you might come across in your very first
# course in calculus

ShowSolution(2*x^(1/2)/Pi^(1/2))

 

Error, (in Student:-Calculus1:-ShowSolution) input expression does not have any incomplete calculus operations

 

NULL

Download fdiff.mw

the simple-minded answer is that the Student[Calculus1] package only contains a limited range of solution methods. Probably sufficient for the sort of problem which might appear in a student textbook, but certainly not representative of the range of problems which Maple can actually handle.

A generally simple method  of determining whether the Student[Calculus1] package *might* be able to come up with step-by-step solution is to use its Hint() command. If this comes up with something "sensible", then use the Rule() command with whatever the Hint() command returns. Keep repeating this approach until a solution is obtained.

The first "toy" example in the attached shows how repeated application of Hint(), Rule() eventually leads to a result.

However the second example in the attached, with your original integrand, shows that the Student[Calculus1] package cannot come up with an approach for the problem. However, assuming that we just don't like the fact that Maple gets a solution anyway, we can apply "wetware" to determine how Maple *might* have achieved this solution.

See the attached (and please read the comments).

  restart:
  with(Student[Calculus1]):
#
# Let us consider a relatively "simple" problem
#
  expr1:= Int( exp(-t)*cos(t), t);
#
# whose solution can be computed immediately
#
  value(expr1);
#
# Assuming we have no idea how this might
# have been computed, let us ask Maple for a
# "hint" about how the original expression
# *might* be transformed, and then apply
# this provided "hint"
#
  h1:= Hint(expr1);
  expr2:= Rule[h1[1]](expr1);
#
# Well this is still not a "solution", so
# let us ask Maple for a "hint" about this
# expression *might* be transformed, and
# then apply this provided "hint"
#
  h2:= Hint(expr2);
  expr3:= Rule[h2[1]](expr2);
#
# Well this is still not a "solution", so
# let us ask Maple for a "hint" about this
# expression *might* be transformed, and
# then apply this provided "hint"
#
  h3:= Hint(expr3);
  expr4:= Rule[h3](expr3);
# Well this is still not a "solution", so
# let us ask Maple for a "hint" about this
# expression *might* be transformed, and
# then apply this provided "hint"
#
  h4:= Hint(expr4);
  expr4:= Rule[h4](expr4);

Int(exp(-t)*cos(t), t)

 

-(1/2)*exp(-t)*cos(t)+(1/2)*exp(-t)*sin(t)

 

[parts, cos(t), -exp(-t), -sin(t), exp(-t), t], [rewrite, exp(-t)*cos(t) = cos(t)/exp(t)]

 

CALCULUS1OBJECT([1, [], []], {t}) = -exp(-t)*cos(t)-(Int(exp(-t)*sin(t), t))

 

[parts, sin(t), -exp(-t), cos(t), exp(-t), t], [rewrite, exp(-t)*sin(t) = sin(t)/exp(t)]

 

CALCULUS1OBJECT([1, [4], []], {t}) = -exp(-t)*cos(t)+exp(-t)*sin(t)+Int(-exp(-t)*cos(t), t)

 

[constantmultiple]

 

CALCULUS1OBJECT([1, [4, 4], []], {t}) = -exp(-t)*cos(t)+exp(-t)*sin(t)-(Int(exp(-t)*cos(t), t))

 

[solve]

 

CALCULUS1OBJECT([1, [], []], {t}) = -(1/2)*exp(-t)*cos(t)+(1/2)*exp(-t)*sin(t)

(1)

#
# So where does applying the above approach get us
# for the OP's original problem, whihc is trivial
# to solve
#
  expr1:= Int( exp(t^2)*cos(t)^3, t);
  value(expr1);
#
# Answer - not very far, because Hint() provides
# nothing.
#
  h1:= Hint(expr1);
#
# So Hint() has failed - what to do now? Well, if we
# really don't want to just accept Maple's solution,
# then we can apply a simple *thought* process.
#
# And my immediate thought is to convert the cos()
# function to an exponential representation, using
# the simple identity that
#
#       cos(x)=(1/2)*exp(I*x)+exp(-I*x))
#
# as in
#
  expr2:= algsubs( cos(t)=(1/2)*(exp(I*t)+exp(-I*t)),
                   Integrand(expr1)
                 )[];
#
# Then expand it into a sum of terms
#
  expr3:=combine( expand(%), power);
#
# Now split the sum of terms into a list
#
  expr4:=[op(expr3)];
#
# Integrate each of these terms individually
#
  expr5:=int~(expr4, t);
#
# Then just add them
#
  add(expr5);

Int(exp(t^2)*cos(t)^3, t)

 

-((1/16)*I)*Pi^(1/2)*exp(9/4)*erf(I*t+3/2)-((3/16)*I)*Pi^(1/2)*exp(1/4)*erf(I*t+1/2)-((3/16)*I)*Pi^(1/2)*exp(1/4)*erf(I*t-1/2)-((1/16)*I)*Pi^(1/2)*exp(9/4)*erf(I*t-3/2)

 

[]

 

(1/8)*exp(t^2)*(exp(I*t)+exp(-I*t))^3

 

(1/8)*exp(t^2+(3*I)*t)+(3/8)*exp(t^2+I*t)+(3/8)*exp(t^2-I*t)+(1/8)*exp(t^2-(3*I)*t)

 

[(1/8)*exp(t^2+(3*I)*t), (3/8)*exp(t^2+I*t), (3/8)*exp(t^2-I*t), (1/8)*exp(t^2-(3*I)*t)]

 

[-((1/16)*I)*Pi^(1/2)*exp(9/4)*erf(I*t-3/2), -((3/16)*I)*Pi^(1/2)*exp(1/4)*erf(I*t-1/2), -((3/16)*I)*Pi^(1/2)*exp(1/4)*erf(I*t+1/2), -((1/16)*I)*Pi^(1/2)*exp(9/4)*erf(I*t+3/2)]

 

-((1/16)*I)*Pi^(1/2)*exp(9/4)*erf(I*t+3/2)-((3/16)*I)*Pi^(1/2)*exp(1/4)*erf(I*t+1/2)-((3/16)*I)*Pi^(1/2)*exp(1/4)*erf(I*t-1/2)-((1/16)*I)*Pi^(1/2)*exp(9/4)*erf(I*t-3/2)

(2)

 

Download doInt.mw

 

demonstrate that neither of  -10.00000000, -8.309536192 can be a solution, so why is solve() producing them?

May be a "numerical issue", but (as far as I can tell) changing the setting of 'Digits' doesn't help :-(

See the attached

  restart;
  Digits:=50:
  expr:=(20 + 20*T + 2*T*(T + 1))*exp(-T) - 10*exp(-2*T) - 2*T - 10. = 0:
#
# Obviously "wrong" solution - no idea why!
#
  ans:=[solve(expr)];
  seq( eval(expr, T=j), j in ans);
#
# Correct(?) solution
#
  plot(lhs(expr), T=-1..10);
  Student[Calculus1]:-Roots(expr, T);

[-10.000000000000000000000000000000000000000000000000, -8.3095361920742886203363304165196635491159255530497, 0.]

 

-4851651944.0979027796910683054154055868463898894485 = 0, -165133630.04379737471898021342327666927866953761632 = 0, 0. = 0

 

 

[0., 1.4114548230186136235259521276376942930399327684067]

(1)

 

 

Download rts.mw

as in

  restart;
  w:=y-x/(c-x):
  p1:=plots:-implicitplot(eval(w,c=2),x=-6..9,y=-3..3, signchange=false):
  p2:=plots:-implicitplot(eval(w,c=2),x=-6..10,y=-3..3, signchange=false):
  p3:=plots:-implicitplot(eval(w,c=2),x=-6..11,y=-3..3, signchange=false):
  [p1,p2,p3]

which produces

 

From the help page for the ':-' operator, (emphasis added);

The binary member selection operator  ":-", is used to name module exports outside the scope in which they are defined. It is a non-commutative, left-associative operator. The first (left) operand is evaluated and must evaluate to a module. The second (right) operand is not evaluated and must be a symbol exported by the module to which the first operand evaluates.

So, for example, the following will "work"

restart;
ode := diff(y(x),x)=a0+a1*y(x)+a2*y(x)^2+a3*y(x)^3;
name_of_package:=DETools;
name_of_package:-abelsol(ode);;

and apply elementwise, as in the "toy" example contained in the attached.

  M:=<1,"", 2;"", 2, 3;4,"","">;
  f:=x-> `if`(x="", 0, x):
  G:=f~(M);

 

Matrix(3, 3, {(1, 1) = 1, (1, 2) = "", (1, 3) = 2, (2, 1) = "", (2, 2) = 2, (2, 3) = 3, (3, 1) = 4, (3, 2) = "", (3, 3) = ""})

 

f := proc (x) options operator, arrow; `if`(x = "", 0, x) end proc

 

Matrix(%id = 36893488148490754396)

(1)

 

Download toyMat.mw

 

about whether you want the k-th power of  a matrix, or the k-th power of the elements of a matrix. For a diagonal matrix, these amount to the same thing. For a non-diagonal matrix, they are different.

See the attached

  restart;
  with(LinearAlgebra):
  DD := <<6, 0> | <0, -1>>:
  DD^~k;
  MatrixPower(DD,k);
  BB := <<6, 1> | <2, -1>>:
  MatrixPower(BB,k);
  BB^~k

Matrix(2, 2, {(1, 1) = 6^k, (1, 2) = 0, (2, 1) = 0, (2, 2) = (-1)^k})

 

Matrix(2, 2, {(1, 1) = 6^k, (1, 2) = 0, (2, 1) = 0, (2, 2) = (-1)^k})

 

Matrix(2, 2, {(1, 1) = (-(7/114)*sqrt(57)+1/2)*(5/2-(1/2)*sqrt(57))^k+(1/114)*(5/2+(1/2)*sqrt(57))^k*(7*sqrt(57)+57), (1, 2) = -(-(2/57)*(5/2+(1/2)*sqrt(57))^k+(2/57)*(5/2-(1/2)*sqrt(57))^k)*sqrt(57), (2, 1) = -(-(1/57)*(5/2+(1/2)*sqrt(57))^k+(1/57)*(5/2-(1/2)*sqrt(57))^k)*sqrt(57), (2, 2) = ((7/114)*sqrt(57)+1/2)*(5/2-(1/2)*sqrt(57))^k+(-(7/114)*sqrt(57)+1/2)*(5/2+(1/2)*sqrt(57))^k})

 

Matrix(%id = 36893488147929398916)

(1)

 


 

Download matPow.mw

for the CartesianProduct command very carefully - at one point it states

CartesianProduct accepts a sequence of graphs as its arguments and returns the Cartesian product of those graphs.  If V1 is the set of vertices of G1, and V2 the set of vertices of G2, then the set of vertices of the Cartesian product G of G1 and G2, is the set V1 X V2.  If u is a vertex in G1, and v a vertex in G2, then the vertices in G are labeled by u:v.  For u1,u2 in V1 and v1,v2 in V2, the edge (u1:v1, u2:v2) is in G iff u1 is adjacent to u2 and v1 = v2 or u1 = u2 and v1 is adjacent to v2.

Isn't that what you want?

 

Everything in Maple is case-sensitive so there is no command COEFF(). There is however a command coeff()), and this appears to work for all (rational) powers.

See the "toy" example in the attached

  restart;
  expr:=add( n^3*x^n, n=-5..5);
  seq( coeff(expr, x, j ), j=-5..5);

-125/x^5-64/x^4-27/x^3-8/x^2-1/x+x+8*x^2+27*x^3+64*x^4+125*x^5

 

-125, -64, -27, -8, -1, 0, 1, 8, 27, 64, 125

(1)

 

Download coeffs.mw

A couple of observations

  1. the first equation 'eq' does not depend on the solution variables you specify. It is a simple cubic with three solutions
  2. the other three equations are symmetric with respect to permutation of the variables x1, x2,  x3. There are six such permutations, so six solution should be obtained
  3. The combination of the above means that 18 solutions are possible

See the attached

  restart:
  eq:= a*x^3 + b*x^2 + c*x + d;
  s:= x1 + x2 + x3 = -b/a;
  sp:= x1*x2 + x1*x3 + x2*x3 = c/a;
  p:= x1*x2*x3 = -d/a;
  [ solve
    ( eval
      ( [eq, p, s, sp],
        [a=1, b=-5, c=6, d=-1]
      ),
      {x, x1, x2, x3},
      explicit
    )
  ]:
  evalf~(%):
  map( fnormal,  %, 1e-08):
  map( simplify, %, zero);
  numelems(%);

a*x^3+b*x^2+c*x+d

 

x1+x2+x3 = -b/a

 

x1*x2+x1*x3+x2*x3 = c/a

 

x1*x2*x3 = -d/a

 

[{x = 3.246979605, x1 = .198062252, x2 = 3.246979605, x3 = 1.554958144}, {x = .1980622645, x1 = 1.554958131, x2 = .1980622645, x3 = 3.246979604}, {x = 1.554958132, x1 = 3.246979604, x2 = 1.554958132, x3 = .198062263}, {x = 3.246979605, x1 = 1.554958144, x2 = 3.246979605, x3 = .198062252}, {x = .1980622645, x1 = 3.246979604, x2 = .1980622645, x3 = 1.554958131}, {x = 1.554958132, x1 = .198062263, x2 = 1.554958132, x3 = 3.246979604}, {x = 3.246979605, x1 = .198062252, x2 = 1.554958144, x3 = 3.246979605}, {x = .1980622645, x1 = 1.554958131, x2 = 3.246979604, x3 = .1980622645}, {x = 1.554958132, x1 = 3.246979604, x2 = .198062263, x3 = 1.554958132}, {x = 3.246979605, x1 = 3.246979605, x2 = 1.554958144, x3 = .198062252}, {x = .1980622645, x1 = .1980622645, x2 = 3.246979604, x3 = 1.554958131}, {x = 1.554958132, x1 = 1.554958132, x2 = .198062263, x3 = 3.246979604}, {x = 3.246979605, x1 = 1.554958144, x2 = .198062252, x3 = 3.246979605}, {x = .1980622645, x1 = 3.246979604, x2 = 1.554958131, x3 = .1980622645}, {x = 1.554958132, x1 = .198062263, x2 = 3.246979604, x3 = 1.554958132}, {x = 3.246979605, x1 = 3.246979605, x2 = .198062252, x3 = 1.554958144}, {x = .1980622645, x1 = .1980622645, x2 = 1.554958131, x3 = 3.246979604}, {x = 1.554958132, x1 = 1.554958132, x2 = 3.246979604, x3 = .198062263}]

 

18

(1)

 


 

Download oddSol.mw

with local/global variables.

For more explanation and solution see the attached

  restart:
  EQ := proc(x, t)
             local y, m, func;
             m := 3;
             func := m*x^2 + t*y^2 = 4;
        end proc:
  eq := EQ(x, 10);
#
# The variable 'y' is local to the procedure EQ
# and thus cannot be assigned a value outside that
# procedure
#
# So both of the following will "fail"
#
  eval(eq, [x=1, y=1]);
  plots:-implicitplot(eq, x= -1 .. 1, y=0..2);
#
# One can "force" force conversion to a global
# variable although this would generally be
# considered bad programming practice
#
  eq:=convert(eq, `global`);
#
# But now both of the following will succeed
#
  eval(eq, [x=1, y=1]);
  plots:-implicitplot(eq, x= -1 .. 1, y=0..2);
 

3*x^2+10*y^2 = 4

 

10*y^2+3 = 4

 

Error, (in plots/implicitplot) invalid input: the following extra unknowns were found in the input expression: {y}

 

3*x^2+10*y^2 = 4

 

13 = 4

 

 

#
# One *could* declare the variable 'y' to be global
# within the  procedure. Although again, declaring
# variables as "global" may be considered bad
# programming practice
#
  restart:
  EQ := proc(x, t)
             local m, func;
             global y;
             m := 3;
             func := m*x^2 + t*y^2 = 4;
        end proc:
  eq := EQ(x, 10);
  eval(eq, [x=1, y=1]);
  plots:-implicitplot(eq, x= -1 .. 1, y=0..2);

3*x^2+10*y^2 = 4

 

13 = 4

 

 

#
# How I would probably do it, avoiding issues
# with local/global
#
  restart:
  EQ := proc(x, y, t)
             local m, func;
             m := 3;
             func := m*x^2 + t*y^2 = 4;
        end proc:
  plots:-implicitplot(EQ(x, y, 10), x= -1 .. 1,  y=0..2);

 

 

Download locGlob.mw

 

 

 - can't really see why you would want to do this

A couple of alternatives are shown in the attached

  restart;
  e:= a/b:
  f:= unapply(unapply(e, [indets(e, name)[]]), x);
  f(z)(A,B);
  g:= unapply(unapply(e, [indets(e, name)[]]));
  g()(A,B);
  h:= unapply(e, [indets(e, name)[]]);
  h(A,B);
  

proc (x) options operator, arrow; proc (a, b) options operator, arrow; a/b end proc end proc

 

A/B

 

proc () options operator, arrow; proc (a, b) options operator, arrow; a/b end proc end proc

 

A/B

 

proc (a, b) options operator, arrow; a/b end proc

 

A/B

(1)

 

Download oddProc2.mw

You have written

C[1, 1] = 0.0998238989835086492681507032141,

I assume you mean

C[1, 1](0) = 0.0998238989835086492681507032141

and the same for alll other initial conditions, in which case I get the attached which executes correctly

restart; with(plots)

sysM := [diff(C[1, 1](t), t) = -(3/16)*Pi*(C[2, 1](t)+4*C[1, 1](t)), diff(C[1, 2](t), t) = -5*Pi*(C[2, 2](t)+4*C[1, 2](t)), diff(C[1, 3](t), t) = -(945/4)*Pi*(C[2, 3](t)+4*C[1, 3](t)), diff(C[2, 1](t), t) = -(1/8)*Pi*(C[3, 1](t)+6*C[2, 1](t)+6*C[1, 1](t)), diff(C[2, 2](t), t) = -10.4719755119659774615421446110*C[3, 2](t)-62.8318530717958647692528676658*C[2, 2](t)-62.8318530717958647692528676658*C[1, 2](t)-2.38361014507273884349657421134*10^15*ZETA[1](t)*C[1, 1](t), diff(C[2, 3](t), t) = -494.800842940392435057866332869*C[3, 3](t)-2968.80505764235461034719799721*C[2, 3](t)-2968.80505764235461034719799721*C[1, 3](t)-1.35954060126371030332767566128*10^16*ZETA[1](t)*C[1, 2](t)-1.35954060126371030332767566128*10^16*ZETA[2](t)*C[1, 1](t), diff(C[3, 1](t), t) = -(3/8)*Pi*(2*C[3, 1](t)+3*C[2, 1](t)), diff(C[3, 2](t), t) = -62.8318530717958647692528676658*C[3, 2](t)-94.2477796076937971538793014986*C[2, 2](t)-1.12625579354686910355213131486*10^17*ZETA[1](t)*C[2, 1](t), diff(C[3, 3](t), t) = -2968.80505764235461034719799721*C[3, 3](t)-4453.20758646353191552079699581*C[2, 3](t)-6.42382934097103118322326749959*10^17*ZETA[1](t)*C[2, 2](t)-6.42382934097103118322326749959*10^17*ZETA[2](t)*C[2, 1](t), diff(ZETA[1](t), t) = -(1/3)*C[2, 1](t), diff(ZETA[2](t), t) = -(1/3)*C[2, 2](t), diff(ZETA[3](t), t) = -(1/3)*C[2, 3](t)]

ICS := [C[1, 1](0) = 0.998238989835086492681507032141e-1, C[1, 2](0) = -0.137051161872492529218951625903e-1, C[1, 3](0) = -0.629146365720807620696267926206e-2, C[2, 1](0) = 0.923300332435106257640735267282e-1, C[2, 2](0) = -0.126762613568515069966837491839e-1, C[2, 3](0) = -0.581915808273854734025727975244e-2, C[3, 1](0) = -0.190143920352772604950256237747e-1, C[3, 2](0) = 0.261054171122321128306306984717e-2, C[3, 3](0) = 0.119839394846335068333097530793e-2, ZETA[1](0) = .464598743230343884242076682299, ZETA[2](0) = .429720916976380440572769663279, ZETA[3](0) = -0.884964696113332752036741498040e-1]

sol := dsolve([sysM[], ICS[]], numeric)

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 .. 28, [( 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) = 12, (2) = 12, (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, (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.2353547231998853e-18, (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..12, {(1) = 0.998238989835086e-1, (2) = -0.137051161872493e-1, (3) = -0.629146365720808e-2, (4) = 0.923300332435106e-1, (5) = -0.126762613568515e-1, (6) = -0.581915808273855e-2, (7) = -0.190143920352773e-1, (8) = 0.261054171122321e-2, (9) = 0.119839394846335e-2, (10) = .464598743230344, (11) = .42972091697638, (12) = -0.884964696113333e-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..12, {(1) = .1, (2) = .1, (3) = .1, (4) = .1, (5) = .1, (6) = .1, (7) = .1, (8) = .1, (9) = .1, (10) = .1, (11) = .1, (12) = .1}, datatype = float[8], order = C_order), Array(1..12, {(1) = .0, (2) = .0, (3) = .0, (4) = .0, (5) = .0, (6) = .0, (7) = .0, (8) = .0, (9) = .0, (10) = .0, (11) = .0, (12) = .0}, datatype = float[8], order = C_order), Array(1..12, {(1) = .0, (2) = .0, (3) = .0, (4) = .0, (5) = .0, (6) = .0, (7) = .0, (8) = .0, (9) = .0, (10) = .0, (11) = .0, (12) = .0}, datatype = float[8], order = C_order), Array(1..12, {(1) = .0, (2) = .0, (3) = .0, (4) = .0, (5) = .0, (6) = .0, (7) = .0, (8) = .0, (9) = .0, (10) = .0, (11) = .0, (12) = .0}, datatype = float[8], order = C_order), Array(1..12, {(1) = .0, (2) = .0, (3) = .0, (4) = .0, (5) = .0, (6) = .0, (7) = .0, (8) = .0, (9) = .0, (10) = .0, (11) = .0, (12) = .0}, datatype = float[8], order = C_order), Array(1..12, 1..12, {(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, (1, 10) = .0, (1, 11) = .0, (1, 12) = .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, (2, 10) = .0, (2, 11) = .0, (2, 12) = .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, (3, 10) = .0, (3, 11) = .0, (3, 12) = .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, (4, 10) = .0, (4, 11) = .0, (4, 12) = .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, (5, 10) = .0, (5, 11) = .0, (5, 12) = .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, (6, 10) = .0, (6, 11) = .0, (6, 12) = .0, (7, 1) = .0, (7, 2) = .0, (7, 3) = .0, (7, 4) = .0, (7, 5) = .0, (7, 6) = .0, (7, 7) = .0, (7, 8) = .0, (7, 9) = .0, (7, 10) = .0, (7, 11) = .0, (7, 12) = .0, (8, 1) = .0, (8, 2) = .0, (8, 3) = .0, (8, 4) = .0, (8, 5) = .0, (8, 6) = .0, (8, 7) = .0, (8, 8) = .0, (8, 9) = .0, (8, 10) = .0, (8, 11) = .0, (8, 12) = .0, (9, 1) = .0, (9, 2) = .0, (9, 3) = .0, (9, 4) = .0, (9, 5) = .0, (9, 6) = .0, (9, 7) = .0, (9, 8) = .0, (9, 9) = .0, (9, 10) = .0, (9, 11) = .0, (9, 12) = .0, (10, 1) = .0, (10, 2) = .0, (10, 3) = .0, (10, 4) = .0, (10, 5) = .0, (10, 6) = .0, (10, 7) = .0, (10, 8) = .0, (10, 9) = .0, (10, 10) = .0, (10, 11) = .0, (10, 12) = .0, (11, 1) = .0, (11, 2) = .0, (11, 3) = .0, (11, 4) = .0, (11, 5) = .0, (11, 6) = .0, (11, 7) = .0, (11, 8) = .0, (11, 9) = .0, (11, 10) = .0, (11, 11) = .0, (11, 12) = .0, (12, 1) = .0, (12, 2) = .0, (12, 3) = .0, (12, 4) = .0, (12, 5) = .0, (12, 6) = .0, (12, 7) = .0, (12, 8) = .0, (12, 9) = .0, (12, 10) = .0, (12, 11) = .0, (12, 12) = .0}, datatype = float[8], order = C_order), Array(1..12, 1..12, {(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, (1, 10) = .0, (1, 11) = .0, (1, 12) = .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, (2, 10) = .0, (2, 11) = .0, (2, 12) = .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, (3, 10) = .0, (3, 11) = .0, (3, 12) = .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, (4, 10) = .0, (4, 11) = .0, (4, 12) = .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, (5, 10) = .0, (5, 11) = .0, (5, 12) = .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, (6, 10) = .0, (6, 11) = .0, (6, 12) = .0, (7, 1) = .0, (7, 2) = .0, (7, 3) = .0, (7, 4) = .0, (7, 5) = .0, (7, 6) = .0, (7, 7) = .0, (7, 8) = .0, (7, 9) = .0, (7, 10) = .0, (7, 11) = .0, (7, 12) = .0, (8, 1) = .0, (8, 2) = .0, (8, 3) = .0, (8, 4) = .0, (8, 5) = .0, (8, 6) = .0, (8, 7) = .0, (8, 8) = .0, (8, 9) = .0, (8, 10) = .0, (8, 11) = .0, (8, 12) = .0, (9, 1) = .0, (9, 2) = .0, (9, 3) = .0, (9, 4) = .0, (9, 5) = .0, (9, 6) = .0, (9, 7) = .0, (9, 8) = .0, (9, 9) = .0, (9, 10) = .0, (9, 11) = .0, (9, 12) = .0, (10, 1) = .0, (10, 2) = .0, (10, 3) = .0, (10, 4) = .0, (10, 5) = .0, (10, 6) = .0, (10, 7) = .0, (10, 8) = .0, (10, 9) = .0, (10, 10) = .0, (10, 11) = .0, (10, 12) = .0, (11, 1) = .0, (11, 2) = .0, (11, 3) = .0, (11, 4) = .0, (11, 5) = .0, (11, 6) = .0, (11, 7) = .0, (11, 8) = .0, (11, 9) = .0, (11, 10) = .0, (11, 11) = .0, (11, 12) = .0, (12, 1) = .0, (12, 2) = .0, (12, 3) = .0, (12, 4) = .0, (12, 5) = .0, (12, 6) = .0, (12, 7) = .0, (12, 8) = .0, (12, 9) = .0, (12, 10) = .0, (12, 11) = .0, (12, 12) = .0}, datatype = float[8], order = C_order), Array(1..12, {(1) = .0, (2) = .0, (3) = .0, (4) = .0, (5) = .0, (6) = .0, (7) = .0, (8) = .0, (9) = .0, (10) = .0, (11) = .0, (12) = .0}, datatype = float[8], order = C_order), Array(1..12, 1..12, {(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, (1, 10) = .0, (1, 11) = .0, (1, 12) = .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, (2, 10) = .0, (2, 11) = .0, (2, 12) = .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, (3, 10) = .0, (3, 11) = .0, (3, 12) = .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, (4, 10) = .0, (4, 11) = .0, (4, 12) = .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, (5, 10) = .0, (5, 11) = .0, (5, 12) = .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, (6, 10) = .0, (6, 11) = .0, (6, 12) = .0, (7, 1) = .0, (7, 2) = .0, (7, 3) = .0, (7, 4) = .0, (7, 5) = .0, (7, 6) = .0, (7, 7) = .0, (7, 8) = .0, (7, 9) = .0, (7, 10) = .0, (7, 11) = .0, (7, 12) = .0, (8, 1) = .0, (8, 2) = .0, (8, 3) = .0, (8, 4) = .0, (8, 5) = .0, (8, 6) = .0, (8, 7) = .0, (8, 8) = .0, (8, 9) = .0, (8, 10) = .0, (8, 11) = .0, (8, 12) = .0, (9, 1) = .0, (9, 2) = .0, (9, 3) = .0, (9, 4) = .0, (9, 5) = .0, (9, 6) = .0, (9, 7) = .0, (9, 8) = .0, (9, 9) = .0, (9, 10) = .0, (9, 11) = .0, (9, 12) = .0, (10, 1) = .0, (10, 2) = .0, (10, 3) = .0, (10, 4) = .0, (10, 5) = .0, (10, 6) = .0, (10, 7) = .0, (10, 8) = .0, (10, 9) = .0, (10, 10) = .0, (10, 11) = .0, (10, 12) = .0, (11, 1) = .0, (11, 2) = .0, (11, 3) = .0, (11, 4) = .0, (11, 5) = .0, (11, 6) = .0, (11, 7) = .0, (11, 8) = .0, (11, 9) = .0, (11, 10) = .0, (11, 11) = .0, (11, 12) = .0, (12, 1) = .0, (12, 2) = .0, (12, 3) = .0, (12, 4) = .0, (12, 5) = .0, (12, 6) = .0, (12, 7) = .0, (12, 8) = .0, (12, 9) = .0, (12, 10) = .0, (12, 11) = .0, (12, 12) = .0}, datatype = float[8], order = C_order), Array(1..12, 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, (10, 1) = .0, (10, 2) = .0, (10, 3) = .0, (10, 4) = .0, (10, 5) = .0, (10, 6) = .0, (11, 1) = .0, (11, 2) = .0, (11, 3) = .0, (11, 4) = .0, (11, 5) = .0, (11, 6) = .0, (12, 1) = .0, (12, 2) = .0, (12, 3) = .0, (12, 4) = .0, (12, 5) = .0, (12, 6) = .0}, datatype = float[8], order = C_order), Array(1..12, {(1) = 0, (2) = 0, (3) = 0, (4) = 0, (5) = 0, (6) = 0, (7) = 0, (8) = 0, (9) = 0, (10) = 0, (11) = 0, (12) = 0}, datatype = integer[8]), Array(1..12, {(1) = .0, (2) = .0, (3) = .0, (4) = .0, (5) = .0, (6) = .0, (7) = .0, (8) = .0, (9) = .0, (10) = .0, (11) = .0, (12) = .0}, datatype = float[8], order = C_order), Array(1..12, {(1) = .0, (2) = .0, (3) = .0, (4) = .0, (5) = .0, (6) = .0, (7) = .0, (8) = .0, (9) = .0, (10) = .0, (11) = .0, (12) = .0}, datatype = float[8], order = C_order), Array(1..12, {(1) = .0, (2) = .0, (3) = .0, (4) = .0, (5) = .0, (6) = .0, (7) = .0, (8) = .0, (9) = .0, (10) = .0, (11) = .0, (12) = .0}, datatype = float[8], order = C_order), Array(1..12, {(1) = .0, (2) = .0, (3) = .0, (4) = .0, (5) = .0, (6) = .0, (7) = .0, (8) = .0, (9) = .0, (10) = .0, (11) = .0, (12) = .0}, datatype = float[8], order = C_order), Array(1..12, {(1) = .0, (2) = .0, (3) = .0, (4) = .0, (5) = .0, (6) = .0, (7) = .0, (8) = .0, (9) = .0, (10) = .0, (11) = .0, (12) = .0}, datatype = float[8], order = C_order), Array(1..24, {(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, (19) = .0, (20) = .0, (21) = .0, (22) = .0, (23) = .0, (24) = .0}, datatype = float[8], order = C_order), Array(1..12, {(1) = 0, (2) = 0, (3) = 0, (4) = 0, (5) = 0, (6) = 0, (7) = 0, (8) = 0, (9) = 0, (10) = 0, (11) = 0, (12) = 0}, datatype = integer[8])]), ( 8 ) = ([Array(1..12, {(1) = 0.998238989835086e-1, (2) = -0.137051161872493e-1, (3) = -0.629146365720808e-2, (4) = 0.923300332435106e-1, (5) = -0.126762613568515e-1, (6) = -0.581915808273855e-2, (7) = -0.190143920352773e-1, (8) = 0.261054171122321e-2, (9) = 0.119839394846335e-2, (10) = .464598743230344, (11) = .42972091697638, (12) = -0.884964696113333e-1}, datatype = float[8], order = C_order), Array(1..12, {(1) = .0, (2) = .0, (3) = .0, (4) = .0, (5) = .0, (6) = .0, (7) = .0, (8) = .0, (9) = .0, (10) = .0, (11) = .0, (12) = .0}, datatype = float[8], order = C_order), Array(1..12, {(1) = -.2895913996763696, (2) = 1.0602360943774856, (3) = 22.997115612306096, (4) = -.44528510209078204, (5) = -110547209583016.8, (6) = -496627098018920.7, (7) = -.28151966766357545, (8) = -4831233864687141., (9) = -21704045380777476., (10) = -0.30776677747836868e-1, (11) = 0.42254204522838325e-2, (12) = 0.19397193609128499e-2}, datatype = float[8], order = C_order), 0, 0]), ( 11 ) = (Array(1..6, 0..12, {(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, (1, 10) = .0, (1, 11) = .0, (1, 12) = .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, (2, 10) = .0, (2, 11) = .0, (2, 12) = .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, (3, 10) = .0, (3, 11) = .0, (3, 12) = .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, (4, 10) = .0, (4, 11) = .0, (4, 12) = .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, (5, 10) = .0, (5, 11) = .0, (5, 12) = .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, (6, 10) = .0, (6, 11) = .0, (6, 12) = .0}, datatype = float[8], order = C_order)), ( 10 ) = ([proc (N, X, Y, YP) option `[Y[1] = C[1,1](t), Y[2] = C[1,2](t), Y[3] = C[1,3](t), Y[4] = C[2,1](t), Y[5] = C[2,2](t), Y[6] = C[2,3](t), Y[7] = C[3,1](t), Y[8] = C[3,2](t), Y[9] = C[3,3](t), Y[10] = ZETA[1](t), Y[11] = ZETA[2](t), Y[12] = ZETA[3](t)]`; YP[1] := -.589048622548086*Y[4]-2.35619449019235*Y[1]; YP[2] := -15.7079632679490*Y[5]-62.8318530717960*Y[2]; YP[3] := -742.201264410588*Y[6]-2968.80505764236*Y[3]; YP[4] := -.392699081698724*Y[7]-2.35619449019234*Y[4]-2.35619449019234*Y[1]; YP[5] := -10.4719755119659774615421446110*Y[8]-62.8318530717958647692528676658*Y[5]-62.8318530717958647692528676658*Y[2]-0.2383610145e16*Y[10]*Y[1]; YP[6] := -494.800842940392435057866332869*Y[9]-2968.80505764235461034719799721*Y[6]-2968.80505764235461034719799721*Y[3]-0.1359540601e17*Y[10]*Y[2]-0.1359540601e17*Y[11]*Y[1]; YP[7] := -2.35619449019234*Y[7]-3.53429173528851*Y[4]; YP[8] := -62.8318530717958647692528676658*Y[8]-94.2477796076937971538793014986*Y[5]-0.1126255794e18*Y[10]*Y[4]; YP[9] := -2968.80505764235461034719799721*Y[9]-4453.20758646353191552079699581*Y[6]-0.6423829341e18*Y[10]*Y[5]-0.6423829341e18*Y[11]*Y[4]; YP[10] := -(1/3)*Y[4]; YP[11] := -(1/3)*Y[5]; YP[12] := -(1/3)*Y[6]; 0 end proc, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]), ( 13 ) = (), ( 12 ) = (), ( 15 ) = ("rkf45"), ( 14 ) = ([0, 0]), ( 18 ) = ([]), ( 19 ) = (0), ( 16 ) = ([0, 0, 0, 0, 0, []]), ( 17 ) = ([proc (N, X, Y, YP) option `[Y[1] = C[1,1](t), Y[2] = C[1,2](t), Y[3] = C[1,3](t), Y[4] = C[2,1](t), Y[5] = C[2,2](t), Y[6] = C[2,3](t), Y[7] = C[3,1](t), Y[8] = C[3,2](t), Y[9] = C[3,3](t), Y[10] = ZETA[1](t), Y[11] = ZETA[2](t), Y[12] = ZETA[3](t)]`; YP[1] := -.589048622548086*Y[4]-2.35619449019235*Y[1]; YP[2] := -15.7079632679490*Y[5]-62.8318530717960*Y[2]; YP[3] := -742.201264410588*Y[6]-2968.80505764236*Y[3]; YP[4] := -.392699081698724*Y[7]-2.35619449019234*Y[4]-2.35619449019234*Y[1]; YP[5] := -10.4719755119659774615421446110*Y[8]-62.8318530717958647692528676658*Y[5]-62.8318530717958647692528676658*Y[2]-0.2383610145e16*Y[10]*Y[1]; YP[6] := -494.800842940392435057866332869*Y[9]-2968.80505764235461034719799721*Y[6]-2968.80505764235461034719799721*Y[3]-0.1359540601e17*Y[10]*Y[2]-0.1359540601e17*Y[11]*Y[1]; YP[7] := -2.35619449019234*Y[7]-3.53429173528851*Y[4]; YP[8] := -62.8318530717958647692528676658*Y[8]-94.2477796076937971538793014986*Y[5]-0.1126255794e18*Y[10]*Y[4]; YP[9] := -2968.80505764235461034719799721*Y[9]-4453.20758646353191552079699581*Y[6]-0.6423829341e18*Y[10]*Y[5]-0.6423829341e18*Y[11]*Y[4]; YP[10] := -(1/3)*Y[4]; YP[11] := -(1/3)*Y[5]; YP[12] := -(1/3)*Y[6]; 0 end proc, -1, 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), ( 28 ) = (0)  ] ))  ] ); _y0 := Array(0..12, {(1) = 0., (2) = 0.998238989835086e-1, (3) = -0.137051161872493e-1, (4) = -0.629146365720808e-2, (5) = 0.923300332435106e-1, (6) = -0.126762613568515e-1, (7) = -0.581915808273855e-2, (8) = -0.190143920352773e-1, (9) = 0.261054171122321e-2, (10) = 0.119839394846335e-2, (11) = .464598743230344, (12) = .429720916976380}); _vmap := array( 1 .. 12, [( 1 ) = (1), ( 2 ) = (2), ( 3 ) = (3), ( 4 ) = (4), ( 5 ) = (5), ( 6 ) = (6), ( 7 ) = (7), ( 9 ) = (9), ( 8 ) = (8), ( 11 ) = (11), ( 10 ) = (10), ( 12 ) = (12)  ] ); _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); _i := false; if _par <> [] then _i := `dsolve/numeric/process_parameters`(_n, _pars, _par, _y0) end if; if _ini <> [] then _i := `dsolve/numeric/process_initial`(_n-_ne, _ini, _y0, _pars, _vmap) or _i end if; if _i then `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 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 elif type(_xin, `=`) and lhs(_xin) = "setdatacallback" then if not type(rhs(_xin), 'nonegint') then error "data callback must be a nonnegative integer (address)" end if; _dtbl[1][28] := rhs(_xin) 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 _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][_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]; 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, _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, C[1, 1](t), C[1, 2](t), C[1, 3](t), C[2, 1](t), C[2, 2](t), C[2, 3](t), C[3, 1](t), C[3, 2](t), C[3, 3](t), ZETA[1](t), ZETA[2](t), ZETA[3](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 := 1; _ndsol := _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)

odeplot(sol, [[t, C[1, 1](t)], [t, C[1, 2](t)], [t, C[1, 3](t)], [t, C[2, 1](t)], [t, C[2, 2](t)], [t, C[2, 3](t)], [t, C[3, 1](t)], [t, C[3, 2](t)], [t, C[3, 3](t)], [t, ZETA[1](t)], [t, ZETA[2](t)], [t, ZETA[3](t)]], t = 0 .. 2)

 

 

Download bigODE.mw

I just downloaded the latest Physics Updates for Maple 2022, and your problem is solved - see the atached

NB I am aware  that you specified Maple 2021: I'm guesing that the lalest Physics updates *may* not work with that release

Attempt to solve with Heaviside

restart;
interface(version);
Physics:-Version();

`Standard Worksheet Interface, Maple 2022.1, Windows 7, May 26 2022 Build ID 1619613`

 

`The "Physics Updates" version in the MapleCloud is 1257 and is the same as the version installed in this computer, created 2022, June 22, 16:27 hours Pacific Time.`

(1)

de := diff(u(x),x$2) = Heaviside(x - a)*u(x);

diff(diff(u(x), x), x) = Heaviside(x-a)*u(x)

(2)

dsolve fails:

dsolve(de);

u(x) = DESol({diff(diff(_Y(x), x), x)-Heaviside(x-a)*_Y(x)}, {_Y(x)})

(3)

Attempt to solve with piecewise

restart;

de := diff(u(x),x$2) = piecewise(x < a, 0, 1)*u(x);

de := diff(u(x), x, x) = piecewise(x < a, 0, 1)*u(x)

(4)

dsolve(de);

u(x) = piecewise(x < a, _C1*x+_C2, a <= x, ((1/2)*(a-1)*exp(a-x)+(1/2)*(a+1)*exp(x-a))*_C1+((1/2)*exp(a-x)+(1/2)*exp(x-a))*_C2)

(5)

dsolve(de, u(x));

u(x) = piecewise(x < a, _C1*x+_C2, a <= x, ((1/2)*(a-1)*exp(a-x)+(1/2)*(a+1)*exp(x-a))*_C1+((1/2)*exp(a-x)+(1/2)*exp(x-a))*_C2)

(6)

 

 

The solution is easy to calculate by hand

We just solve the (quite trivial) DE over the intervals x < a and x > 0

separately, and patch the two solutions by requiring the continuity

of u(x) and diff(u(x), x) at x = a.  We get

sol := piecewise(x < a,
        x*c[1] + c[2],
        ((a*c[1] + c[1] + c[2])*exp(x))/(2*exp(a)) + ((a*c[1] - c[1] + c[2])*exp(-x))/(2*exp(-a)));

sol := piecewise(x < a, x*c[1]+c[2], (a*c[1]+c[1]+c[2])*exp(x)/(2*exp(a))+(a*c[1]-c[1]+c[2])*exp(-x)/(2*exp(-a)))

(7)

 

Download badHeavi.mw

I answered this question about 6 hours ago, after which both the question and my answer disappeared.

I feel disinclined to give the same lengthy explanations about the required corrections again, so the attached is posted without comment.


 

  restart;
  with(plots):

  eq1:= diff(f(x), x$3)
        +
        1/2*(1-phi)^2.5*(1-phi+phi*rho__s/rho__f1)*(eta*cos__omega + f(x)*sin__omega)*diff(f(x), x$2)
        +
        (1-phi)^2.5*M*sin(alpha)^2*(1-diff(f(x),x))
        +
        (1-phi)^2.5*(1-phi+phi*rho__beta__s/rho__beta__f1)*lambda__T*theta(x);

  eq2:= K__nf/K__f*diff(theta(x),x$2)
        +
        Pr/2*(eta*cos__omega + f(x)*sin__omega)*diff(theta(x),x);
#
# The OP's original list of boundary conditions contained f(1)=0.
# In a 1-dimensional BVP only two boundaries are allowed, so the
# condition f(1)=0 has been "replaced" by the condition
#
#    D(f)(0) = val
#
# A "shooting" method will be used later to determine a value for
# "val", whihc ensures that the original BC ( ie f(1)=0 ) is
# satisfied
#
  bcs := f(0) = 0, D(f)(0) = val, f(10) = 1, theta(0) = 1, theta(10) = 0;

  params:= [ K__f=0.613, K__nf=0.6842, M=1, Pr=6.2,
             alpha=0, eta=0, phi=0.05, rho__f1=997.1,
             rho__s=5200, cos__omega=1, lambda__T=0,
             rho__beta__f1=20939.1, rho__beta__s=997.1,
             sin__omega=1
           ];
#
# So the ODE system with all parameter values inserted is
#
  odesys:=eval([eq1, eq2], params);

diff(diff(diff(f(x), x), x), x)+(1/2)*(1-phi)^2.5*(1-phi+phi*rho__s/rho__f1)*(eta*cos__omega+f(x)*sin__omega)*(diff(diff(f(x), x), x))+(1-phi)^2.5*M*sin(alpha)^2*(1-(diff(f(x), x)))+(1-phi)^2.5*(1-phi+phi*rho__beta__s/rho__beta__f1)*lambda__T*theta(x)

 

K__nf*(diff(diff(theta(x), x), x))/K__f+(1/2)*Pr*(eta*cos__omega+f(x)*sin__omega)*(diff(theta(x), x))

 

f(0) = 0, (D(f))(0) = val, f(10) = 1, theta(0) = 1, theta(10) = 0

 

[K__f = .613, K__nf = .6842, M = 1, Pr = 6.2, alpha = 0, eta = 0, phi = 0.5e-1, rho__f1 = 997.1, rho__s = 5200, cos__omega = 1, lambda__T = 0, rho__beta__f1 = 20939.1, rho__beta__s = 997.1, sin__omega = 1]

 

[diff(diff(diff(f(x), x), x), x)+.5325197465*f(x)*(diff(diff(f(x), x), x)), 1.116150082*(diff(diff(theta(x), x), x))+3.100000000*f(x)*(diff(theta(x), x))]

(1)

#
# The shooting method
#
# This determines a value for D(f)(0) such that the
# OP's original desired condition f(1)=0 is satisfied
#
  SM:= proc( shoot )
             local ans:
             if type( shoot, numeric)
             then ans:= dsolve
                        ( eval
                          ( [odesys[], bcs],
                            val=shoot
                          ),
                          numeric
                        ):
                  return rhs(ans(1)[2]);
             else return 'procname(_passed)'
             fi;
       end proc:
  sv:= fsolve( SM );

-0.1290531852e-1

(2)

#
# Solve the system
#
  sol:= dsolve
        (  eval
           ( [odesys[], bcs],
             val=sv
           ),
          numeric
        ):
#
# Check that when x=1, f(x)=0, which was OP's desired condition
#
  sol(1)[1..2];
#
# Plot the results
#
  odeplot
  ( sol,
    [ [x, f(x)],
      [x, theta(x)]
    ],
    x = 0 .. 10,
    color=[red, blue],
    legend=[typeset(f(x)), typeset(theta(x))]
  );
  

[x = 1., f(x) = HFloat(1.9274077505456594e-12)]

 

 

 


 

Download odeProb.mw

5 6 7 8 9 10 11 Last Page 7 of 191