acer

32353 Reputation

29 Badges

19 years, 331 days
Ontario, Canada

Social Networks and Content at Maplesoft.com

MaplePrimes Activity


These are replies submitted by acer

@Gibberish123 I see. So for your (n,e) RSA scheme you'll need to generate a number num from the plaintext string such that gcd(num,n)=1 , and 0<=num<n .

Once you can do that, the encryption/decryption might be as follows:

restart:

n := 17;
                                       n := 17

L := numtheory[lambda](n);
                                       L := 16

e := 5;
                                       e := 5

evalb( e < L ); # a requirement 
                                        true

evalb( gcd(e, L) = 1); # a requirement
                                        true

d := 1/e mod L;
                                       d := 13

num := 10;  # a number generated from the plaintext string
                                       num := 10

c := num^e mod n;  # the encrypted number
                                        c := 6

c^d mod n;  # decrypting back to num
                                          10

In Maple, for a large n, it is more efficient to use the following command for modular exponentiation (which avoids explicit computation of num^e which happens when executing num^e mod n ).

modp(Power(num,e), n); # efficient modular exponentiation
                                           6

modp(Power(c,d), n); # efficient modular exponentiation  
                                           10

So given a plaintext ascii string you have to use an n, as well as figure out a way to convert from string to number, and pad, so that (at least) the following conditions hold:

evalb( gcd(num,n)=1 and num<n );
                                         true

@tomleslie It looks to me like your suggestion [edited to remove suggestion to divide by N] is the composite trapezoidal rule (or the average of the left and right Riemann sums, but not a mid-point Riemann sum).

It's not a wrong approach. But quite often one might attain a more accurate result with a higher order interpolation scheme (ie. higher than linear). Of course "quite often" relates to the smoothness of well-behavedness of the underlying function approximated by the data.

The OP has not mentioned accuracy or performance, unfortunately.

 

@Gibberish123 Correct, I did not implement RSA encryption/decryption for you. Are you hoping for someone to implement a secure padding scheme such as this? By the way, is this course work?

@Markiyan Hirnyk Your previous comments to your own answer above had,

int(c[2]*(diff(c[1], v)), v = 0 .. 9);

instead of,

int(c[2]*(diff(c[1], v)), v = 1 .. 6);

which you now use. My reply to you is that your initial upper end-point range for the parameter (ie, 9), while integrating, was invalid, which I would have thought was obvious in my previous response. You have not shown how you got "6" here, either.

Let's hope that the OP's data represents a mathematical funcion (one y value only for each x value), and not a more general curve.

@Markiyan Hirnyk I posted an Answer to the OP's question, not a Reply to your Answer. The OP did not supply data, as far as I have noticed, and as I was not replying to your Answer I see no special reason to use the same example data as you have.

And I have also used different interpolating commands from you, yes. I don't think that BSplineCurve is necessarily the right thing to use: there is a reason why the 3rd entry in the list it returns may not be a range (for the parameter denoting x(v) and y(v) ) that extends the curve as a mathematical function y(x) over all the independent data. A B-spline is a beast of another color, and one should be extra careful if trying to force it to cover the full range of the (supplied) independent data if a function (mathematical sense) is desired.

By "independent range as returned by BSplineCurve" of course I mean a range for the parameter (not the independent data variable).


 

restart;

with(CurveFitting):

xydata := Matrix([[0,0],[1,1],[4,9],[6,10],[8,5],[8,3]]):

c := BSplineCurve(xydata,v):

c[3];   # independent range returned by BSplineCurve

v = 3 .. 6

min(xydata[..,1]) .. max(xydata[..,1]); # independent range in xydata

0 .. 8

min(xydata[..,2]) .. max(xydata[..,2]); # dependent range in xydata

0 .. 10

plot(c, view=[min(xydata[..,1]) .. max(xydata[..,1]),
              min(xydata[..,2]) .. max(xydata[..,2])]);

plot([c[1],c[2],v = min(xydata[..,1]) .. max(xydata[..,1])]);

 


 

Download bsplinething.mw

@AmusingYeti I find `tools/genglobal` pretty useful, when I want to generate an unassigned global name. (It gets passed the base of the name.) Sometimes I use PiecewiseTools:-Flatten. But usefulness is in the eye of the beholder.

@sand15 In Maple 2015.2 I had to put the interface calls in separate Execution Groups in order to have it work, using a Worksheet. (That seems fixed in Maple 2016.0, btw.)


 

restart;

kernelopts(version);

`Maple 2015.2, X86 64 WINDOWS, Nov 13 2015, Build ID 1087698`

F:=dsolve({diff(x(t),t)=sin(t),x(0)=-1},{x(t)},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 .. 24, [( 1 ) = (datatype = float[8], order = C_order, storage = rectangular), ( 2 ) = (datatype = float[8], order = C_order, storage = rectangular), ( 3 ) = ([0, 0, 0, Array(1..0, 1..2, {}, datatype = float[8], order = C_order)]), ( 4 ) = (Array(1..54, {(1) = 1, (2) = 1, (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) = 21, (19) = 30000, (20) = 0, (21) = 0, (22) = 1, (23) = 4, (24) = 0, (25) = 1, (26) = 15, (27) = 1, (28) = 0, (29) = 1, (30) = 3, (31) = 3, (32) = 0, (33) = 1, (34) = 0, (35) = 0, (36) = 0, (37) = 0, (38) = 0, (39) = 0, (40) = 0, (41) = 0, (42) = 0, (43) = 1, (44) = 0, (45) = 0, (46) = 0, (47) = 0, (48) = 0, (49) = 0, (50) = 50, (51) = 1, (52) = 0, (53) = 0, (54) = 0}, datatype = integer[8])), ( 5 ) = (Array(1..28, {(1) = .0, (2) = 0.10e-5, (3) = .0, (4) = 0.500001e-14, (5) = .0, (6) = .16695651424556227, (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..1, {(1) = -1.0}, 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..1, {(1) = .1}, datatype = float[8], order = C_order), Array(1..1, {(1) = .0}, datatype = float[8], order = C_order), Array(1..1, {(1) = .0}, datatype = float[8], order = C_order), Array(1..1, {(1) = .0}, datatype = float[8], order = C_order), Array(1..1, {(1) = .0}, datatype = float[8], order = C_order), Array(1..1, 1..1, {(1, 1) = .0}, datatype = float[8], order = C_order), Array(1..1, 1..1, {(1, 1) = .0}, datatype = float[8], order = C_order), Array(1..1, 1..6, {(1, 1) = .0, (1, 2) = .0, (1, 3) = .0, (1, 4) = .0, (1, 5) = .0, (1, 6) = .0}, datatype = float[8], order = C_order), Array(1..1, {(1) = 0}, datatype = integer[8]), Array(1..1, {(1) = .0}, datatype = float[8], order = C_order), Array(1..1, {(1) = .0}, datatype = float[8], order = C_order), Array(1..1, {(1) = .0}, datatype = float[8], order = C_order), Array(1..1, {(1) = .0}, datatype = float[8], order = C_order), Array(1..1, {(1) = .0}, datatype = float[8], order = C_order)]), ( 8 ) = ([Array(1..1, {(1) = -1.0}, datatype = float[8], order = C_order), Array(1..1, {(1) = .0}, datatype = float[8], order = C_order), Array(1..1, {(1) = .0}, datatype = float[8], order = C_order), 0, 0]), ( 11 ) = (Array(1..6, 0..1, {(1, 1) = .0, (2, 0) = .0, (2, 1) = .0, (3, 0) = .0, (3, 1) = .0, (4, 0) = .0, (4, 1) = .0, (5, 0) = .0, (5, 1) = .0, (6, 0) = .0, (6, 1) = .0}, datatype = float[8], order = C_order)), ( 10 ) = ([proc (N, X, Y, YP) option `[Y[1] = x(t)]`; YP[1] := sin(X); 0 end proc, -1, 0, 0, 0, 0, 0, 0]), ( 13 ) = (), ( 12 ) = (), ( 15 ) = ("rkf45"), ( 14 ) = ([0, 0]), ( 18 ) = ([]), ( 19 ) = (0), ( 16 ) = ([0, 0, 0, []]), ( 17 ) = ([proc (N, X, Y, YP) option `[Y[1] = x(t)]`; YP[1] := sin(X); 0 end proc, -1, 0, 0, 0, 0, 0, 0]), ( 22 ) = (0), ( 23 ) = (0), ( 20 ) = ([]), ( 21 ) = (0), ( 24 ) = (0)  ] ))  ] ); _y0 := Array(0..1, {(1) = 0.}); _vmap := array( 1 .. 1, [( 1 ) = (1)  ] ); _x0 := _dtbl[1][5][5]; _n := _dtbl[1][4][1]; _ne := _dtbl[1][4][3]; _nd := _dtbl[1][4][4]; _nv := _dtbl[1][4][16]; if not type(_xout, 'numeric') then if member(_xout, ["start", "left", "right"]) then if _Env_smart_dsolve_numeric = true or _dtbl[1][4][10] = 1 then if _xout = "left" then if type(_dtbl[2], 'table') then return _dtbl[2][5][1] end if elif _xout = "right" then if type(_dtbl[3], 'table') then return _dtbl[3][5][1] end if end if end if; return _dtbl[1][5][5] elif _xout = "method" then return _dtbl[1][15] elif _xout = "storage" then return evalb(_dtbl[1][4][10] = 1) elif _xout = "leftdata" then if not type(_dtbl[2], 'array') then return NULL else return eval(_dtbl[2]) end if elif _xout = "rightdata" then if not type(_dtbl[3], 'array') then return NULL else return eval(_dtbl[3]) end if elif _xout = "enginedata" then return eval(_dtbl[1]) elif _xout = "enginereset" then _dtbl[2] := evaln(_dtbl[2]); _dtbl[3] := evaln(_dtbl[3]); return NULL elif _xout = "initial" then return procname(_y0[0]) elif _xout = "laxtol" then return _dtbl[`if`(member(_dtbl[4], {2, 3}), _dtbl[4], 1)][5][18] elif _xout = "numfun" then return `if`(member(_dtbl[4], {2, 3}), _dtbl[_dtbl[4]][4][18], 0) elif _xout = "parameters" then return [seq(_y0[_n+_i], _i = 1 .. nops(_pars))] elif _xout = "initial_and_parameters" then return procname(_y0[0]), [seq(_y0[_n+_i], _i = 1 .. nops(_pars))] elif _xout = "last" then if _dtbl[4] <> 2 and _dtbl[4] <> 3 or _x0-_dtbl[_dtbl[4]][5][1] = 0. then error "no information is available on last computed point" else _xout := _dtbl[_dtbl[4]][5][1] end if elif _xout = "function" then if _dtbl[1][4][33]-2. = 0 then return eval(_dtbl[1][10], 1) else return eval(_dtbl[1][10][1], 1) end if elif _xout = "map" then return copy(_vmap) elif type(_xin, `=`) and type(rhs(_xin), 'list') and member(lhs(_xin), {"initial", "parameters", "initial_and_parameters"}) then _ini, _par := [], []; if lhs(_xin) = "initial" then _ini := rhs(_xin) elif lhs(_xin) = "parameters" then _par := rhs(_xin) elif select(type, rhs(_xin), `=`) <> [] then _par, _ini := selectremove(type, rhs(_xin), `=`) elif nops(rhs(_xin)) < nops(_pars)+1 then error "insufficient data for specification of initial and parameters" else _par := rhs(_xin)[-nops(_pars) .. -1]; _ini := rhs(_xin)[1 .. -nops(_pars)-1] end if; _xout := lhs(_xout); if _par <> [] then `dsolve/numeric/process_parameters`(_n, _pars, _par, _y0) end if; if _ini <> [] then `dsolve/numeric/process_initial`(_n-_ne, _ini, _y0, _pars, _vmap) end if; `dsolve/numeric/SC/reinitialize`(_dtbl, _y0, _n, procname, _pars); if _Env_smart_dsolve_numeric = true and type(_y0[0], 'numeric') and _dtbl[1][4][10] <> 1 then procname("right") := _y0[0]; procname("left") := _y0[0] end if; if _xout = "initial" then return [_y0[0], seq(_y0[_vmap[_i]], _i = 1 .. _n-_ne)] elif _xout = "parameters" then return [seq(_y0[_n+_i], _i = 1 .. nops(_pars))] else return [_y0[0], seq(_y0[_vmap[_i]], _i = 1 .. _n-_ne)], [seq(_y0[_n+_i], _i = 1 .. nops(_pars))] end if elif _xin = "eventstop" then if _nv = 0 then error "this solution has no events" end if; _i := _dtbl[4]; if _i <> 2 and _i <> 3 then return 0 end if; if _dtbl[_i][4][10] = 1 and assigned(_dtbl[5-_i]) and _dtbl[_i][4][9] < 100 and 100 <= _dtbl[5-_i][4][9] then _i := 5-_i; _dtbl[4] := _i; _j := round(_dtbl[_i][4][17]); return round(_dtbl[_i][3][1][_j, 1]) elif 100 <= _dtbl[_i][4][9] then _j := round(_dtbl[_i][4][17]); return round(_dtbl[_i][3][1][_j, 1]) else return 0 end if elif _xin = "eventstatus" then if _nv = 0 then error "this solution has no events" end if; _i := [selectremove(proc (a) options operator, arrow; _dtbl[1][3][1][a, 7] = 1 end proc, {seq(_j, _j = 1 .. round(_dtbl[1][3][1][_nv+1, 1]))})]; return ':-enabled' = _i[1], ':-disabled' = _i[2] elif _xin = "eventclear" then if _nv = 0 then error "this solution has no events" end if; _i := _dtbl[4]; if _i <> 2 and _i <> 3 then error "no events to clear" end if; if _dtbl[_i][4][10] = 1 and assigned(_dtbl[5-_i]) and _dtbl[_i][4][9] < 100 and 100 < _dtbl[5-_i][4][9] then _dtbl[4] := 5-_i; _i := 5-_i end if; if _dtbl[_i][4][9] < 100 then error "no events to clear" elif _nv < _dtbl[_i][4][9]-100 then error "event error condition cannot be cleared" else _j := _dtbl[_i][4][9]-100; if irem(round(_dtbl[_i][3][1][_j, 4]), 2) = 1 then error "retriggerable events cannot be cleared" end if; _j := round(_dtbl[_i][3][1][_j, 1]); for _k to _nv do if _dtbl[_i][3][1][_k, 1] = _j then if _dtbl[_i][3][1][_k, 2] = 3 then error "range events cannot be cleared" end if; _dtbl[_i][3][1][_k, 8] := _dtbl[_i][3][1][_nv+1, 8] end if end do; _dtbl[_i][4][17] := 0; _dtbl[_i][4][9] := 0; if _dtbl[1][4][10] = 1 then if _i = 2 then try procname(procname("left")) catch:  end try else try procname(procname("right")) catch:  end try end if end if end if; return  elif type(_xin, `=`) and member(lhs(_xin), {"eventdisable", "eventenable"}) then if _nv = 0 then error "this solution has no events" end if; if type(rhs(_xin), {('list')('posint'), ('set')('posint')}) then _i := {op(rhs(_xin))} elif type(rhs(_xin), 'posint') then _i := {rhs(_xin)} else error "event identifiers must be integers in the range 1..%1", round(_dtbl[1][3][1][_nv+1, 1]) end if; if select(proc (a) options operator, arrow; _nv < a end proc, _i) <> {} then error "event identifiers must be integers in the range 1..%1", round(_dtbl[1][3][1][_nv+1, 1]) end if; _k := {}; for _j to _nv do if member(round(_dtbl[1][3][1][_j, 1]), _i) then _k := `union`(_k, {_j}) end if end do; _i := _k; if lhs(_xin) = "eventdisable" then _dtbl[4] := 0; _j := [evalb(assigned(_dtbl[2]) and member(_dtbl[2][4][17], _i)), evalb(assigned(_dtbl[3]) and member(_dtbl[3][4][17], _i))]; for _k in _i do _dtbl[1][3][1][_k, 7] := 0; if assigned(_dtbl[2]) then _dtbl[2][3][1][_k, 7] := 0 end if; if assigned(_dtbl[3]) then _dtbl[3][3][1][_k, 7] := 0 end if end do; if _j[1] then for _k to _nv+1 do if _k <= _nv and not type(_dtbl[2][3][4][_k, 1], 'undefined') then userinfo(3, {'events', 'eventreset'}, `reinit #2, event code `, _k, ` to defined init `, _dtbl[2][3][4][_k, 1]); _dtbl[2][3][1][_k, 8] := _dtbl[2][3][4][_k, 1] elif _dtbl[2][3][1][_k, 2] = 0 and irem(iquo(round(_dtbl[2][3][1][_k, 4]), 32), 2) = 1 then userinfo(3, {'events', 'eventreset'}, `reinit #2, event code `, _k, ` to rate hysteresis init `, _dtbl[2][5][24]); _dtbl[2][3][1][_k, 8] := _dtbl[2][5][24] elif _dtbl[2][3][1][_k, 2] = 0 and irem(iquo(round(_dtbl[2][3][1][_k, 4]), 2), 2) = 0 then userinfo(3, {'events', 'eventreset'}, `reinit #2, event code `, _k, ` to initial init `, _x0); _dtbl[2][3][1][_k, 8] := _x0 else userinfo(3, {'events', 'eventreset'}, `reinit #2, event code `, _k, ` to fireinitial init `, _x0-1); _dtbl[2][3][1][_k, 8] := _x0-1 end if end do; _dtbl[2][4][17] := 0; _dtbl[2][4][9] := 0; if _dtbl[1][4][10] = 1 then procname(procname("left")) end if end if; if _j[2] then for _k to _nv+1 do if _k <= _nv and not type(_dtbl[3][3][4][_k, 2], 'undefined') then userinfo(3, {'events', 'eventreset'}, `reinit #3, event code `, _k, ` to defined init `, _dtbl[3][3][4][_k, 2]); _dtbl[3][3][1][_k, 8] := _dtbl[3][3][4][_k, 2] elif _dtbl[3][3][1][_k, 2] = 0 and irem(iquo(round(_dtbl[3][3][1][_k, 4]), 32), 2) = 1 then userinfo(3, {'events', 'eventreset'}, `reinit #3, event code `, _k, ` to rate hysteresis init `, _dtbl[3][5][24]); _dtbl[3][3][1][_k, 8] := _dtbl[3][5][24] elif _dtbl[3][3][1][_k, 2] = 0 and irem(iquo(round(_dtbl[3][3][1][_k, 4]), 2), 2) = 0 then userinfo(3, {'events', 'eventreset'}, `reinit #3, event code `, _k, ` to initial init `, _x0); _dtbl[3][3][1][_k, 8] := _x0 else userinfo(3, {'events', 'eventreset'}, `reinit #3, event code `, _k, ` to fireinitial init `, _x0+1); _dtbl[3][3][1][_k, 8] := _x0+1 end if end do; _dtbl[3][4][17] := 0; _dtbl[3][4][9] := 0; if _dtbl[1][4][10] = 1 then procname(procname("right")) end if end if else for _k in _i do _dtbl[1][3][1][_k, 7] := 1 end do; _dtbl[2] := evaln(_dtbl[2]); _dtbl[3] := evaln(_dtbl[3]); _dtbl[4] := 0; if _dtbl[1][4][10] = 1 then if _x0 <= procname("right") then try procname(procname("right")) catch:  end try end if; if procname("left") <= _x0 then try procname(procname("left")) catch:  end try end if end if end if; return  elif type(_xin, `=`) and lhs(_xin) = "eventfired" then if not type(rhs(_xin), 'list') then error "'eventfired' must be specified as a list" end if; if _nv = 0 then error "this solution has no events" end if; if _dtbl[4] <> 2 and _dtbl[4] <> 3 then error "'direction' must be set prior to calling/setting 'eventfired'" end if; _i := _dtbl[4]; _val := NULL; if not assigned(_EnvEventRetriggerWarned) then _EnvEventRetriggerWarned := false end if; for _k in rhs(_xin) do if type(_k, 'integer') then _src := _k elif type(_k, 'integer' = 'anything') and type(evalf(rhs(_k)), 'numeric') then _k := lhs(_k) = evalf[max(Digits, 18)](rhs(_k)); _src := lhs(_k) else error "'eventfired' entry is not valid: %1", _k end if; if _src < 1 or round(_dtbl[1][3][1][_nv+1, 1]) < _src then error "event identifiers must be integers in the range 1..%1", round(_dtbl[1][3][1][_nv+1, 1]) end if; _src := {seq(`if`(_dtbl[1][3][1][_j, 1]-_src = 0., _j, NULL), _j = 1 .. _nv)}; if nops(_src) <> 1 then error "'eventfired' can only be set/queried for root-finding events and time/interval events" end if; _src := _src[1]; if _dtbl[1][3][1][_src, 2] <> 0. and _dtbl[1][3][1][_src, 2]-2. <> 0. then error "'eventfired' can only be set/queried for root-finding events and time/interval events" elif irem(round(_dtbl[1][3][1][_src, 4]), 2) = 1 then if _EnvEventRetriggerWarned = false then WARNING(`'eventfired' has no effect on events that retrigger`) end if; _EnvEventRetriggerWarned := true end if; if _dtbl[_i][3][1][_src, 2] = 0 and irem(iquo(round(_dtbl[_i][3][1][_src, 4]), 32), 2) = 1 then _val := _val, undefined elif type(_dtbl[_i][3][4][_src, _i-1], 'undefined') or _i = 2 and _dtbl[2][3][1][_src, 8] < _dtbl[2][3][4][_src, 1] or _i = 3 and _dtbl[3][3][4][_src, 2] < _dtbl[3][3][1][_src, 8] then _val := _val, _dtbl[_i][3][1][_src, 8] else _val := _val, _dtbl[_i][3][4][_src, _i-1] end if; if type(_k, `=`) then if _dtbl[_i][3][1][_src, 2] = 0 and irem(iquo(round(_dtbl[_i][3][1][_src, 4]), 32), 2) = 1 then error "cannot set event code for a rate hysteresis event" end if; userinfo(3, {'events', 'eventreset'}, `manual set event code `, _src, ` to value `, rhs(_k)); _dtbl[_i][3][1][_src, 8] := rhs(_k); _dtbl[_i][3][4][_src, _i-1] := rhs(_k) end if end do; return [_val] elif type(_xin, `=`) and lhs(_xin) = "direction" then if not member(rhs(_xin), {-1, 1, ':-left', ':-right'}) then error "'direction' must be specified as either '1' or 'right' (positive) or '-1' or 'left' (negative)" end if; _src := `if`(_dtbl[4] = 2, -1, `if`(_dtbl[4] = 3, 1, undefined)); _i := `if`(member(rhs(_xin), {1, ':-right'}), 3, 2); _dtbl[4] := _i; _dtbl[_i] := `dsolve/numeric/SC/IVPdcopy`(_dtbl[1], `if`(assigned(_dtbl[_i]), _dtbl[_i], NULL)); if 0 < _nv then for _j to _nv+1 do if _j <= _nv and not type(_dtbl[_i][3][4][_j, _i-1], 'undefined') then userinfo(3, {'events', 'eventreset'}, `reinit #4, event code `, _j, ` to defined init `, _dtbl[_i][3][4][_j, _i-1]); _dtbl[_i][3][1][_j, 8] := _dtbl[_i][3][4][_j, _i-1] elif _dtbl[_i][3][1][_j, 2] = 0 and irem(iquo(round(_dtbl[_i][3][1][_j, 4]), 32), 2) = 1 then userinfo(3, {'events', 'eventreset'}, `reinit #4, event code `, _j, ` to rate hysteresis init `, _dtbl[_i][5][24]); _dtbl[_i][3][1][_j, 8] := _dtbl[_i][5][24] elif _dtbl[_i][3][1][_j, 2] = 0 and irem(iquo(round(_dtbl[_i][3][1][_j, 4]), 2), 2) = 0 then userinfo(3, {'events', 'eventreset'}, `reinit #4, event code `, _j, ` to initial init `, _x0); _dtbl[_i][3][1][_j, 8] := _x0 else userinfo(3, {'events', 'eventreset'}, `reinit #4, event code `, _j, ` to fireinitial init `, _x0-2*_i+5.0); _dtbl[_i][3][1][_j, 8] := _x0-2*_i+5.0 end if end do end if; return _src elif _xin = "eventcount" then if _dtbl[1][3][1] = 0 or _dtbl[4] <> 2 and _dtbl[4] <> 3 then return 0 else return round(_dtbl[_dtbl[4]][3][1][_nv+1, 12]) end if else return "procname" end if end if; if _xout = _x0 then return [_x0, seq(evalf(_dtbl[1][6][_vmap[_i]]), _i = 1 .. _n-_ne)] end if; _i := `if`(_x0 <= _xout, 3, 2); if _xin = "last" and 0 < _dtbl[_i][4][9] and _dtbl[_i][4][9] < 100 then _dat := eval(_dtbl[_i], 2); _j := _dat[4][20]; return [_dat[11][_j, 0], seq(_dat[11][_j, _vmap[_i]], _i = 1 .. _n-_ne-_nd), seq(_dat[8][1][_vmap[_i]], _i = _n-_ne-_nd+1 .. _n-_ne)] end if; if not type(_dtbl[_i], 'array') then _dtbl[_i] := `dsolve/numeric/SC/IVPdcopy`(_dtbl[1], `if`(assigned(_dtbl[_i]), _dtbl[_i], NULL)); if 0 < _nv then for _j to _nv+1 do if _j <= _nv and not type(_dtbl[_i][3][4][_j, _i-1], 'undefined') then userinfo(3, {'events', 'eventreset'}, `reinit #5, event code `, _j, ` to defined init `, _dtbl[_i][3][4][_j, _i-1]); _dtbl[_i][3][1][_j, 8] := _dtbl[_i][3][4][_j, _i-1] elif _dtbl[_i][3][1][_j, 2] = 0 and irem(iquo(round(_dtbl[_i][3][1][_j, 4]), 32), 2) = 1 then userinfo(3, {'events', 'eventreset'}, `reinit #5, event code `, _j, ` to rate hysteresis init `, _dtbl[_i][5][24]); _dtbl[_i][3][1][_j, 8] := _dtbl[_i][5][24] elif _dtbl[_i][3][1][_j, 2] = 0 and irem(iquo(round(_dtbl[_i][3][1][_j, 4]), 2), 2) = 0 then userinfo(3, {'events', 'eventreset'}, `reinit #5, event code `, _j, ` to initial init `, _x0); _dtbl[_i][3][1][_j, 8] := _x0 else userinfo(3, {'events', 'eventreset'}, `reinit #5, event code `, _j, ` to fireinitial init `, _x0-2*_i+5.0); _dtbl[_i][3][1][_j, 8] := _x0-2*_i+5.0 end if end do end if end if; if _xin <> "last" then if 0 < 0 then if `dsolve/numeric/checkglobals`(op(_dtbl[1][14]), _pars, _n, _y0) then `dsolve/numeric/SC/reinitialize`(_dtbl, _y0, _n, procname, _pars, _i) end if end if; if _dtbl[1][4][7] = 0 then error "parameters must be initialized before solution can be computed" end if end if; _dat := eval(_dtbl[_i], 2); _dtbl[4] := _i; try _src := `dsolve/numeric/SC/IVPrun`(_dat, _xout) catch: userinfo(2, `dsolve/debug`, print(`Exception in solnproc:`, [lastexception][2 .. -1])); error  end try; if _src = 0 and 100 < _dat[4][9] then _val := _dat[3][1][_nv+1, 8] else _val := _dat[11][_dat[4][20], 0] end if; if _src <> 0 or _dat[4][9] <= 0 then _dtbl[1][5][1] := _xout else _dtbl[1][5][1] := _val end if; if _i = 3 and _val < _xout then Rounding := -infinity; if _dat[4][9] = 1 then error "cannot evaluate the solution further right of %1, probably a singularity", evalf[8](_val) elif _dat[4][9] = 2 then error "cannot evaluate the solution further right of %1, maxfun limit exceeded (see ?dsolve,maxfun for details)", evalf[8](_val) elif _dat[4][9] = 3 then if _dat[4][25] = 3 then error "cannot evaluate the solution past the initial point, problem may be initially singular or improperly set up" else error "cannot evaluate the solution past the initial point, problem may be complex, initially singular or improperly set up" end if elif _dat[4][9] = 4 then error "cannot evaluate the solution further right of %1, accuracy goal cannot be achieved with specified 'minstep'", evalf[8](_val) elif _dat[4][9] = 5 then error "cannot evaluate the solution further right of %1, too many step failures, tolerances may be too loose for problem", evalf[8](_val) elif _dat[4][9] = 6 then error "cannot evaluate the solution further right of %1, cannot downgrade delay storage for problems with delay derivative order > 1, try increasing delaypts", evalf[8](_val) elif _dat[4][9] = 10 then error "cannot evaluate the solution further right of %1, interrupt requested", evalf[8](_val) elif 100 < _dat[4][9] then if _dat[4][9]-100 = _nv+1 then error "constraint projection failure on event at t=%1", evalf[8](_val) elif _dat[4][9]-100 = _nv+2 then error "index-1 and derivative evaluation failure on event at t=%1", evalf[8](_val) elif _dat[4][9]-100 = _nv+3 then error "maximum number of event iterations reached (%1) at t=%2", round(_dat[3][1][_nv+1, 3]), evalf[8](_val) else if _Env_dsolve_nowarnstop <> true then `dsolve/numeric/warning`(StringTools:-FormatMessage("cannot evaluate the solution further right of %1, event #%2 triggered a halt", evalf[8](_val), round(_dat[3][1][_dat[4][9]-100, 1]))) end if; Rounding := 'nearest'; _xout := _val end if else error "cannot evaluate the solution further right of %1", evalf[8](_val) end if elif _i = 2 and _xout < _val then Rounding := infinity; if _dat[4][9] = 1 then error "cannot evaluate the solution further left of %1, probably a singularity", evalf[8](_val) elif _dat[4][9] = 2 then error "cannot evaluate the solution further left of %1, maxfun limit exceeded (see ?dsolve,maxfun for details)", evalf[8](_val) elif _dat[4][9] = 3 then if _dat[4][25] = 3 then error "cannot evaluate the solution past the initial point, problem may be initially singular or improperly set up" else error "cannot evaluate the solution past the initial point, problem may be complex, initially singular or improperly set up" end if elif _dat[4][9] = 4 then error "cannot evaluate the solution further left of %1, accuracy goal cannot be achieved with specified 'minstep'", evalf[8](_val) elif _dat[4][9] = 5 then error "cannot evaluate the solution further left of %1, too many step failures, tolerances may be too loose for problem", evalf[8](_val) elif _dat[4][9] = 6 then error "cannot evaluate the solution further left of %1, cannot downgrade delay storage for problems with delay derivative order > 1, try increasing delaypts", evalf[8](_val) elif _dat[4][9] = 10 then error "cannot evaluate the solution further right of %1, interrupt requested", evalf[8](_val) elif 100 < _dat[4][9] then if _dat[4][9]-100 = _nv+1 then error "constraint projection failure on event at t=%1", evalf[8](_val) elif _dat[4][9]-100 = _nv+2 then error "index-1 and derivative evaluation failure on event at t=%1", evalf[8](_val) elif _dat[4][9]-100 = _nv+3 then error "maximum number of event iterations reached (%1) at t=%2", round(_dat[3][1][_nv+1, 3]), evalf[8](_val) else if _Env_dsolve_nowarnstop <> true then `dsolve/numeric/warning`(StringTools:-FormatMessage("cannot evaluate the solution further left of %1, event #%2 triggered a halt", evalf[8](_val), round(_dat[3][1][_dat[4][9]-100, 1]))) end if; Rounding := 'nearest'; _xout := _val end if else error "cannot evaluate the solution further left of %1", evalf[8](_val) end if end if; if _EnvInFsolve = true then _dig := _dat[4][26]; _dat[4][26] := _EnvDSNumericSaveDigits; _Env_dsolve_SC_native := true; if _dat[4][25] = 1 then _i := 1; _dat[4][25] := 2 else _i := _dat[4][25] end if; _val := `dsolve/numeric/SC/IVPval`(_dat, _xout, _src); _dat[4][25] := _i; _dat[4][26] := _dig; [_xout, seq(_val[_vmap[_i]], _i = 1 .. _n-_ne)] else Digits := _dat[4][26]; _val := `dsolve/numeric/SC/IVPval`(eval(_dat, 2), _xout, _src); [_xout, seq(_val[_vmap[_i]], _i = 1 .. _n-_ne)] end if end proc, (2) = Array(0..0, {}), (3) = [t, x(t)], (4) = []}); _vars := _dat[3]; _pars := map(rhs, _dat[4]); _n := nops(_vars)-1; _solnproc := _dat[1]; if not type(_xout, 'numeric') then if member(x_rkf45, ["start", 'start', "method", 'method', "left", 'left', "right", 'right', "leftdata", "rightdata", "enginedata", "eventstop", 'eventstop', "eventclear", 'eventclear', "eventstatus", 'eventstatus', "eventcount", 'eventcount', "laxtol", 'laxtol', "numfun", 'numfun', NULL]) then _res := _solnproc(convert(x_rkf45, 'string')); if 1 < nops([_res]) then return _res elif type(_res, 'array') then return eval(_res, 1) elif _res <> "procname" then return _res end if elif member(x_rkf45, ["last", 'last', "initial", 'initial', "parameters", 'parameters', "initial_and_parameters", 'initial_and_parameters', NULL]) then _xout := convert(x_rkf45, 'string'); _res := _solnproc(_xout); if _xout = "parameters" then return [seq(_pars[_i] = _res[_i], _i = 1 .. nops(_pars))] elif _xout = "initial_and_parameters" then return [seq(_vars[_i+1] = [_res][1][_i+1], _i = 0 .. _n), seq(_pars[_i] = [_res][2][_i], _i = 1 .. nops(_pars))] else return [seq(_vars[_i+1] = _res[_i+1], _i = 0 .. _n)] end if elif type(_xout, `=`) and member(lhs(_xout), ["initial", 'initial', "parameters", 'parameters', "initial_and_parameters", 'initial_and_parameters', NULL]) then _xout := convert(lhs(x_rkf45), 'string') = rhs(x_rkf45); if type(rhs(_xout), 'list') then _res := _solnproc(_xout) else error "initial and/or parameter values must be specified in a list" end if; if lhs(_xout) = "initial" then return [seq(_vars[_i+1] = _res[_i+1], _i = 0 .. _n)] elif lhs(_xout) = "parameters" then return [seq(_pars[_i] = _res[_i], _i = 1 .. nops(_pars))] else return [seq(_vars[_i+1] = [_res][1][_i+1], _i = 0 .. _n), seq(_pars[_i] = [_res][2][_i], _i = 1 .. nops(_pars))] end if elif type(_xout, `=`) and member(lhs(_xout), ["eventdisable", 'eventdisable', "eventenable", 'eventenable', "eventfired", 'eventfired', "direction", 'direction', NULL]) then return _solnproc(convert(lhs(x_rkf45), 'string') = rhs(x_rkf45)) elif _xout = "solnprocedure" then return eval(_solnproc) elif _xout = "sysvars" then return _vars end if; if procname <> unknown then return ('procname')(x_rkf45) else _ndsol; _ndsol := pointto(_dat[2][0]); return ('_ndsol')(x_rkf45) end if end if; try _res := _solnproc(_xout); [seq(_vars[_i+1] = _res[_i+1], _i = 0 .. _n)] catch: error  end try end proc

#plots:-odeplot(F); # -cos(t)

str:="C:\\\\TEMP\\foo.txt":
interface('prettyprint'=1):
interface('verboseproc'=3):

writeto(str);
print(F);
writeto('terminal'):

interface('prettyprint'=3):
interface('verboseproc'=1):

eval(F);

proc (x_rkf45) local _res, _dat, _vars, _solnproc, _xout, _ndsol, _pars, _n, _i; option `Copyright (c) 2000 by Waterloo Maple Inc. All rights reserved.`; if 1 < nargs then error "invalid input: too many arguments" end if; _EnvDSNumericSaveDigits := Digits; Digits := 15; if _EnvInFsolve = true then _xout := evalf[_EnvDSNumericSaveDigits](x_rkf45) else _xout := evalf(x_rkf45) end if; _dat := Array(1..4, {(1) = proc (_xin) local _xout, _dtbl, _dat, _vmap, _x0, _y0, _val, _dig, _n, _ne, _nd, _nv, _pars, _ini, _par, _i, _j, _k, _src; option `Copyright (c) 2002 by Waterloo Maple Inc. All rights reserved.`; table( [( "complex" ) = false ] ) _xout := _xin; _pars := []; _dtbl := array( 1 .. 4, [( 1 ) = (array( 1 .. 24, [( 1 ) = (datatype = float[8], order = C_order, storage = rectangular), ( 2 ) = (datatype = float[8], order = C_order, storage = rectangular), ( 3 ) = ([0, 0, 0, Array(1..0, 1..2, {}, datatype = float[8], order = C_order)]), ( 4 ) = (Array(1..54, {(1) = 1, (2) = 1, (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) = 21, (19) = 30000, (20) = 0, (21) = 0, (22) = 1, (23) = 4, (24) = 0, (25) = 1, (26) = 15, (27) = 1, (28) = 0, (29) = 1, (30) = 3, (31) = 3, (32) = 0, (33) = 1, (34) = 0, (35) = 0, (36) = 0, (37) = 0, (38) = 0, (39) = 0, (40) = 0, (41) = 0, (42) = 0, (43) = 1, (44) = 0, (45) = 0, (46) = 0, (47) = 0, (48) = 0, (49) = 0, (50) = 50, (51) = 1, (52) = 0, (53) = 0, (54) = 0}, datatype = integer[8])), ( 5 ) = (Array(1..28, {(1) = .0, (2) = 0.10e-5, (3) = .0, (4) = 0.500001e-14, (5) = .0, (6) = .16695651424556227, (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..1, {(1) = -1.0}, 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..1, {(1) = .1}, datatype = float[8], order = C_order), Array(1..1, {(1) = .0}, datatype = float[8], order = C_order), Array(1..1, {(1) = .0}, datatype = float[8], order = C_order), Array(1..1, {(1) = .0}, datatype = float[8], order = C_order), Array(1..1, {(1) = .0}, datatype = float[8], order = C_order), Array(1..1, 1..1, {(1, 1) = .0}, datatype = float[8], order = C_order), Array(1..1, 1..1, {(1, 1) = .0}, datatype = float[8], order = C_order), Array(1..1, 1..6, {(1, 1) = .0, (1, 2) = .0, (1, 3) = .0, (1, 4) = .0, (1, 5) = .0, (1, 6) = .0}, datatype = float[8], order = C_order), Array(1..1, {(1) = 0}, datatype = integer[8]), Array(1..1, {(1) = .0}, datatype = float[8], order = C_order), Array(1..1, {(1) = .0}, datatype = float[8], order = C_order), Array(1..1, {(1) = .0}, datatype = float[8], order = C_order), Array(1..1, {(1) = .0}, datatype = float[8], order = C_order), Array(1..1, {(1) = .0}, datatype = float[8], order = C_order)]), ( 8 ) = ([Array(1..1, {(1) = -1.0}, datatype = float[8], order = C_order), Array(1..1, {(1) = .0}, datatype = float[8], order = C_order), Array(1..1, {(1) = .0}, datatype = float[8], order = C_order), 0, 0]), ( 11 ) = (Array(1..6, 0..1, {(1, 1) = .0, (2, 0) = .0, (2, 1) = .0, (3, 0) = .0, (3, 1) = .0, (4, 0) = .0, (4, 1) = .0, (5, 0) = .0, (5, 1) = .0, (6, 0) = .0, (6, 1) = .0}, datatype = float[8], order = C_order)), ( 10 ) = ([proc (N, X, Y, YP) option `[Y[1] = x(t)]`; YP[1] := sin(X); 0 end proc, -1, 0, 0, 0, 0, 0, 0]), ( 13 ) = (), ( 12 ) = (), ( 15 ) = ("rkf45"), ( 14 ) = ([0, 0]), ( 18 ) = ([]), ( 19 ) = (0), ( 16 ) = ([0, 0, 0, []]), ( 17 ) = ([proc (N, X, Y, YP) option `[Y[1] = x(t)]`; YP[1] := sin(X); 0 end proc, -1, 0, 0, 0, 0, 0, 0]), ( 22 ) = (0), ( 23 ) = (0), ( 20 ) = ([]), ( 21 ) = (0), ( 24 ) = (0)  ] ))  ] ); _y0 := Array(0..1, {(1) = 0.}); _vmap := array( 1 .. 1, [( 1 ) = (1)  ] ); _x0 := _dtbl[1][5][5]; _n := _dtbl[1][4][1]; _ne := _dtbl[1][4][3]; _nd := _dtbl[1][4][4]; _nv := _dtbl[1][4][16]; if not type(_xout, 'numeric') then if member(_xout, ["start", "left", "right"]) then if _Env_smart_dsolve_numeric = true or _dtbl[1][4][10] = 1 then if _xout = "left" then if type(_dtbl[2], 'table') then return _dtbl[2][5][1] end if elif _xout = "right" then if type(_dtbl[3], 'table') then return _dtbl[3][5][1] end if end if end if; return _dtbl[1][5][5] elif _xout = "method" then return _dtbl[1][15] elif _xout = "storage" then return evalb(_dtbl[1][4][10] = 1) elif _xout = "leftdata" then if not type(_dtbl[2], 'array') then return NULL else return eval(_dtbl[2]) end if elif _xout = "rightdata" then if not type(_dtbl[3], 'array') then return NULL else return eval(_dtbl[3]) end if elif _xout = "enginedata" then return eval(_dtbl[1]) elif _xout = "enginereset" then _dtbl[2] := evaln(_dtbl[2]); _dtbl[3] := evaln(_dtbl[3]); return NULL elif _xout = "initial" then return procname(_y0[0]) elif _xout = "laxtol" then return _dtbl[`if`(member(_dtbl[4], {2, 3}), _dtbl[4], 1)][5][18] elif _xout = "numfun" then return `if`(member(_dtbl[4], {2, 3}), _dtbl[_dtbl[4]][4][18], 0) elif _xout = "parameters" then return [seq(_y0[_n+_i], _i = 1 .. nops(_pars))] elif _xout = "initial_and_parameters" then return procname(_y0[0]), [seq(_y0[_n+_i], _i = 1 .. nops(_pars))] elif _xout = "last" then if _dtbl[4] <> 2 and _dtbl[4] <> 3 or _x0-_dtbl[_dtbl[4]][5][1] = 0. then error "no information is available on last computed point" else _xout := _dtbl[_dtbl[4]][5][1] end if elif _xout = "function" then if _dtbl[1][4][33]-2. = 0 then return eval(_dtbl[1][10], 1) else return eval(_dtbl[1][10][1], 1) end if elif _xout = "map" then return copy(_vmap) elif type(_xin, `=`) and type(rhs(_xin), 'list') and member(lhs(_xin), {"initial", "parameters", "initial_and_parameters"}) then _ini, _par := [], []; if lhs(_xin) = "initial" then _ini := rhs(_xin) elif lhs(_xin) = "parameters" then _par := rhs(_xin) elif select(type, rhs(_xin), `=`) <> [] then _par, _ini := selectremove(type, rhs(_xin), `=`) elif nops(rhs(_xin)) < nops(_pars)+1 then error "insufficient data for specification of initial and parameters" else _par := rhs(_xin)[-nops(_pars) .. -1]; _ini := rhs(_xin)[1 .. -nops(_pars)-1] end if; _xout := lhs(_xout); if _par <> [] then `dsolve/numeric/process_parameters`(_n, _pars, _par, _y0) end if; if _ini <> [] then `dsolve/numeric/process_initial`(_n-_ne, _ini, _y0, _pars, _vmap) end if; `dsolve/numeric/SC/reinitialize`(_dtbl, _y0, _n, procname, _pars); if _Env_smart_dsolve_numeric = true and type(_y0[0], 'numeric') and _dtbl[1][4][10] <> 1 then procname("right") := _y0[0]; procname("left") := _y0[0] end if; if _xout = "initial" then return [_y0[0], seq(_y0[_vmap[_i]], _i = 1 .. _n-_ne)] elif _xout = "parameters" then return [seq(_y0[_n+_i], _i = 1 .. nops(_pars))] else return [_y0[0], seq(_y0[_vmap[_i]], _i = 1 .. _n-_ne)], [seq(_y0[_n+_i], _i = 1 .. nops(_pars))] end if elif _xin = "eventstop" then if _nv = 0 then error "this solution has no events" end if; _i := _dtbl[4]; if _i <> 2 and _i <> 3 then return 0 end if; if _dtbl[_i][4][10] = 1 and assigned(_dtbl[5-_i]) and _dtbl[_i][4][9] < 100 and 100 <= _dtbl[5-_i][4][9] then _i := 5-_i; _dtbl[4] := _i; _j := round(_dtbl[_i][4][17]); return round(_dtbl[_i][3][1][_j, 1]) elif 100 <= _dtbl[_i][4][9] then _j := round(_dtbl[_i][4][17]); return round(_dtbl[_i][3][1][_j, 1]) else return 0 end if elif _xin = "eventstatus" then if _nv = 0 then error "this solution has no events" end if; _i := [selectremove(proc (a) options operator, arrow; _dtbl[1][3][1][a, 7] = 1 end proc, {seq(_j, _j = 1 .. round(_dtbl[1][3][1][_nv+1, 1]))})]; return ':-enabled' = _i[1], ':-disabled' = _i[2] elif _xin = "eventclear" then if _nv = 0 then error "this solution has no events" end if; _i := _dtbl[4]; if _i <> 2 and _i <> 3 then error "no events to clear" end if; if _dtbl[_i][4][10] = 1 and assigned(_dtbl[5-_i]) and _dtbl[_i][4][9] < 100 and 100 < _dtbl[5-_i][4][9] then _dtbl[4] := 5-_i; _i := 5-_i end if; if _dtbl[_i][4][9] < 100 then error "no events to clear" elif _nv < _dtbl[_i][4][9]-100 then error "event error condition cannot be cleared" else _j := _dtbl[_i][4][9]-100; if irem(round(_dtbl[_i][3][1][_j, 4]), 2) = 1 then error "retriggerable events cannot be cleared" end if; _j := round(_dtbl[_i][3][1][_j, 1]); for _k to _nv do if _dtbl[_i][3][1][_k, 1] = _j then if _dtbl[_i][3][1][_k, 2] = 3 then error "range events cannot be cleared" end if; _dtbl[_i][3][1][_k, 8] := _dtbl[_i][3][1][_nv+1, 8] end if end do; _dtbl[_i][4][17] := 0; _dtbl[_i][4][9] := 0; if _dtbl[1][4][10] = 1 then if _i = 2 then try procname(procname("left")) catch:  end try else try procname(procname("right")) catch:  end try end if end if end if; return  elif type(_xin, `=`) and member(lhs(_xin), {"eventdisable", "eventenable"}) then if _nv = 0 then error "this solution has no events" end if; if type(rhs(_xin), {('list')('posint'), ('set')('posint')}) then _i := {op(rhs(_xin))} elif type(rhs(_xin), 'posint') then _i := {rhs(_xin)} else error "event identifiers must be integers in the range 1..%1", round(_dtbl[1][3][1][_nv+1, 1]) end if; if select(proc (a) options operator, arrow; _nv < a end proc, _i) <> {} then error "event identifiers must be integers in the range 1..%1", round(_dtbl[1][3][1][_nv+1, 1]) end if; _k := {}; for _j to _nv do if member(round(_dtbl[1][3][1][_j, 1]), _i) then _k := `union`(_k, {_j}) end if end do; _i := _k; if lhs(_xin) = "eventdisable" then _dtbl[4] := 0; _j := [evalb(assigned(_dtbl[2]) and member(_dtbl[2][4][17], _i)), evalb(assigned(_dtbl[3]) and member(_dtbl[3][4][17], _i))]; for _k in _i do _dtbl[1][3][1][_k, 7] := 0; if assigned(_dtbl[2]) then _dtbl[2][3][1][_k, 7] := 0 end if; if assigned(_dtbl[3]) then _dtbl[3][3][1][_k, 7] := 0 end if end do; if _j[1] then for _k to _nv+1 do if _k <= _nv and not type(_dtbl[2][3][4][_k, 1], 'undefined') then userinfo(3, {'events', 'eventreset'}, `reinit #2, event code `, _k, ` to defined init `, _dtbl[2][3][4][_k, 1]); _dtbl[2][3][1][_k, 8] := _dtbl[2][3][4][_k, 1] elif _dtbl[2][3][1][_k, 2] = 0 and irem(iquo(round(_dtbl[2][3][1][_k, 4]), 32), 2) = 1 then userinfo(3, {'events', 'eventreset'}, `reinit #2, event code `, _k, ` to rate hysteresis init `, _dtbl[2][5][24]); _dtbl[2][3][1][_k, 8] := _dtbl[2][5][24] elif _dtbl[2][3][1][_k, 2] = 0 and irem(iquo(round(_dtbl[2][3][1][_k, 4]), 2), 2) = 0 then userinfo(3, {'events', 'eventreset'}, `reinit #2, event code `, _k, ` to initial init `, _x0); _dtbl[2][3][1][_k, 8] := _x0 else userinfo(3, {'events', 'eventreset'}, `reinit #2, event code `, _k, ` to fireinitial init `, _x0-1); _dtbl[2][3][1][_k, 8] := _x0-1 end if end do; _dtbl[2][4][17] := 0; _dtbl[2][4][9] := 0; if _dtbl[1][4][10] = 1 then procname(procname("left")) end if end if; if _j[2] then for _k to _nv+1 do if _k <= _nv and not type(_dtbl[3][3][4][_k, 2], 'undefined') then userinfo(3, {'events', 'eventreset'}, `reinit #3, event code `, _k, ` to defined init `, _dtbl[3][3][4][_k, 2]); _dtbl[3][3][1][_k, 8] := _dtbl[3][3][4][_k, 2] elif _dtbl[3][3][1][_k, 2] = 0 and irem(iquo(round(_dtbl[3][3][1][_k, 4]), 32), 2) = 1 then userinfo(3, {'events', 'eventreset'}, `reinit #3, event code `, _k, ` to rate hysteresis init `, _dtbl[3][5][24]); _dtbl[3][3][1][_k, 8] := _dtbl[3][5][24] elif _dtbl[3][3][1][_k, 2] = 0 and irem(iquo(round(_dtbl[3][3][1][_k, 4]), 2), 2) = 0 then userinfo(3, {'events', 'eventreset'}, `reinit #3, event code `, _k, ` to initial init `, _x0); _dtbl[3][3][1][_k, 8] := _x0 else userinfo(3, {'events', 'eventreset'}, `reinit #3, event code `, _k, ` to fireinitial init `, _x0+1); _dtbl[3][3][1][_k, 8] := _x0+1 end if end do; _dtbl[3][4][17] := 0; _dtbl[3][4][9] := 0; if _dtbl[1][4][10] = 1 then procname(procname("right")) end if end if else for _k in _i do _dtbl[1][3][1][_k, 7] := 1 end do; _dtbl[2] := evaln(_dtbl[2]); _dtbl[3] := evaln(_dtbl[3]); _dtbl[4] := 0; if _dtbl[1][4][10] = 1 then if _x0 <= procname("right") then try procname(procname("right")) catch:  end try end if; if procname("left") <= _x0 then try procname(procname("left")) catch:  end try end if end if end if; return  elif type(_xin, `=`) and lhs(_xin) = "eventfired" then if not type(rhs(_xin), 'list') then error "'eventfired' must be specified as a list" end if; if _nv = 0 then error "this solution has no events" end if; if _dtbl[4] <> 2 and _dtbl[4] <> 3 then error "'direction' must be set prior to calling/setting 'eventfired'" end if; _i := _dtbl[4]; _val := NULL; if not assigned(_EnvEventRetriggerWarned) then _EnvEventRetriggerWarned := false end if; for _k in rhs(_xin) do if type(_k, 'integer') then _src := _k elif type(_k, 'integer' = 'anything') and type(evalf(rhs(_k)), 'numeric') then _k := lhs(_k) = evalf[max(Digits, 18)](rhs(_k)); _src := lhs(_k) else error "'eventfired' entry is not valid: %1", _k end if; if _src < 1 or round(_dtbl[1][3][1][_nv+1, 1]) < _src then error "event identifiers must be integers in the range 1..%1", round(_dtbl[1][3][1][_nv+1, 1]) end if; _src := {seq(`if`(_dtbl[1][3][1][_j, 1]-_src = 0., _j, NULL), _j = 1 .. _nv)}; if nops(_src) <> 1 then error "'eventfired' can only be set/queried for root-finding events and time/interval events" end if; _src := _src[1]; if _dtbl[1][3][1][_src, 2] <> 0. and _dtbl[1][3][1][_src, 2]-2. <> 0. then error "'eventfired' can only be set/queried for root-finding events and time/interval events" elif irem(round(_dtbl[1][3][1][_src, 4]), 2) = 1 then if _EnvEventRetriggerWarned = false then WARNING(`'eventfired' has no effect on events that retrigger`) end if; _EnvEventRetriggerWarned := true end if; if _dtbl[_i][3][1][_src, 2] = 0 and irem(iquo(round(_dtbl[_i][3][1][_src, 4]), 32), 2) = 1 then _val := _val, undefined elif type(_dtbl[_i][3][4][_src, _i-1], 'undefined') or _i = 2 and _dtbl[2][3][1][_src, 8] < _dtbl[2][3][4][_src, 1] or _i = 3 and _dtbl[3][3][4][_src, 2] < _dtbl[3][3][1][_src, 8] then _val := _val, _dtbl[_i][3][1][_src, 8] else _val := _val, _dtbl[_i][3][4][_src, _i-1] end if; if type(_k, `=`) then if _dtbl[_i][3][1][_src, 2] = 0 and irem(iquo(round(_dtbl[_i][3][1][_src, 4]), 32), 2) = 1 then error "cannot set event code for a rate hysteresis event" end if; userinfo(3, {'events', 'eventreset'}, `manual set event code `, _src, ` to value `, rhs(_k)); _dtbl[_i][3][1][_src, 8] := rhs(_k); _dtbl[_i][3][4][_src, _i-1] := rhs(_k) end if end do; return [_val] elif type(_xin, `=`) and lhs(_xin) = "direction" then if not member(rhs(_xin), {-1, 1, ':-left', ':-right'}) then error "'direction' must be specified as either '1' or 'right' (positive) or '-1' or 'left' (negative)" end if; _src := `if`(_dtbl[4] = 2, -1, `if`(_dtbl[4] = 3, 1, undefined)); _i := `if`(member(rhs(_xin), {1, ':-right'}), 3, 2); _dtbl[4] := _i; _dtbl[_i] := `dsolve/numeric/SC/IVPdcopy`(_dtbl[1], `if`(assigned(_dtbl[_i]), _dtbl[_i], NULL)); if 0 < _nv then for _j to _nv+1 do if _j <= _nv and not type(_dtbl[_i][3][4][_j, _i-1], 'undefined') then userinfo(3, {'events', 'eventreset'}, `reinit #4, event code `, _j, ` to defined init `, _dtbl[_i][3][4][_j, _i-1]); _dtbl[_i][3][1][_j, 8] := _dtbl[_i][3][4][_j, _i-1] elif _dtbl[_i][3][1][_j, 2] = 0 and irem(iquo(round(_dtbl[_i][3][1][_j, 4]), 32), 2) = 1 then userinfo(3, {'events', 'eventreset'}, `reinit #4, event code `, _j, ` to rate hysteresis init `, _dtbl[_i][5][24]); _dtbl[_i][3][1][_j, 8] := _dtbl[_i][5][24] elif _dtbl[_i][3][1][_j, 2] = 0 and irem(iquo(round(_dtbl[_i][3][1][_j, 4]), 2), 2) = 0 then userinfo(3, {'events', 'eventreset'}, `reinit #4, event code `, _j, ` to initial init `, _x0); _dtbl[_i][3][1][_j, 8] := _x0 else userinfo(3, {'events', 'eventreset'}, `reinit #4, event code `, _j, ` to fireinitial init `, _x0-2*_i+5.0); _dtbl[_i][3][1][_j, 8] := _x0-2*_i+5.0 end if end do end if; return _src elif _xin = "eventcount" then if _dtbl[1][3][1] = 0 or _dtbl[4] <> 2 and _dtbl[4] <> 3 then return 0 else return round(_dtbl[_dtbl[4]][3][1][_nv+1, 12]) end if else return "procname" end if end if; if _xout = _x0 then return [_x0, seq(evalf(_dtbl[1][6][_vmap[_i]]), _i = 1 .. _n-_ne)] end if; _i := `if`(_x0 <= _xout, 3, 2); if _xin = "last" and 0 < _dtbl[_i][4][9] and _dtbl[_i][4][9] < 100 then _dat := eval(_dtbl[_i], 2); _j := _dat[4][20]; return [_dat[11][_j, 0], seq(_dat[11][_j, _vmap[_i]], _i = 1 .. _n-_ne-_nd), seq(_dat[8][1][_vmap[_i]], _i = _n-_ne-_nd+1 .. _n-_ne)] end if; if not type(_dtbl[_i], 'array') then _dtbl[_i] := `dsolve/numeric/SC/IVPdcopy`(_dtbl[1], `if`(assigned(_dtbl[_i]), _dtbl[_i], NULL)); if 0 < _nv then for _j to _nv+1 do if _j <= _nv and not type(_dtbl[_i][3][4][_j, _i-1], 'undefined') then userinfo(3, {'events', 'eventreset'}, `reinit #5, event code `, _j, ` to defined init `, _dtbl[_i][3][4][_j, _i-1]); _dtbl[_i][3][1][_j, 8] := _dtbl[_i][3][4][_j, _i-1] elif _dtbl[_i][3][1][_j, 2] = 0 and irem(iquo(round(_dtbl[_i][3][1][_j, 4]), 32), 2) = 1 then userinfo(3, {'events', 'eventreset'}, `reinit #5, event code `, _j, ` to rate hysteresis init `, _dtbl[_i][5][24]); _dtbl[_i][3][1][_j, 8] := _dtbl[_i][5][24] elif _dtbl[_i][3][1][_j, 2] = 0 and irem(iquo(round(_dtbl[_i][3][1][_j, 4]), 2), 2) = 0 then userinfo(3, {'events', 'eventreset'}, `reinit #5, event code `, _j, ` to initial init `, _x0); _dtbl[_i][3][1][_j, 8] := _x0 else userinfo(3, {'events', 'eventreset'}, `reinit #5, event code `, _j, ` to fireinitial init `, _x0-2*_i+5.0); _dtbl[_i][3][1][_j, 8] := _x0-2*_i+5.0 end if end do end if end if; if _xin <> "last" then if 0 < 0 then if `dsolve/numeric/checkglobals`(op(_dtbl[1][14]), _pars, _n, _y0) then `dsolve/numeric/SC/reinitialize`(_dtbl, _y0, _n, procname, _pars, _i) end if end if; if _dtbl[1][4][7] = 0 then error "parameters must be initialized before solution can be computed" end if end if; _dat := eval(_dtbl[_i], 2); _dtbl[4] := _i; try _src := `dsolve/numeric/SC/IVPrun`(_dat, _xout) catch: userinfo(2, `dsolve/debug`, print(`Exception in solnproc:`, [lastexception][2 .. -1])); error  end try; if _src = 0 and 100 < _dat[4][9] then _val := _dat[3][1][_nv+1, 8] else _val := _dat[11][_dat[4][20], 0] end if; if _src <> 0 or _dat[4][9] <= 0 then _dtbl[1][5][1] := _xout else _dtbl[1][5][1] := _val end if; if _i = 3 and _val < _xout then Rounding := -infinity; if _dat[4][9] = 1 then error "cannot evaluate the solution further right of %1, probably a singularity", evalf[8](_val) elif _dat[4][9] = 2 then error "cannot evaluate the solution further right of %1, maxfun limit exceeded (see ?dsolve,maxfun for details)", evalf[8](_val) elif _dat[4][9] = 3 then if _dat[4][25] = 3 then error "cannot evaluate the solution past the initial point, problem may be initially singular or improperly set up" else error "cannot evaluate the solution past the initial point, problem may be complex, initially singular or improperly set up" end if elif _dat[4][9] = 4 then error "cannot evaluate the solution further right of %1, accuracy goal cannot be achieved with specified 'minstep'", evalf[8](_val) elif _dat[4][9] = 5 then error "cannot evaluate the solution further right of %1, too many step failures, tolerances may be too loose for problem", evalf[8](_val) elif _dat[4][9] = 6 then error "cannot evaluate the solution further right of %1, cannot downgrade delay storage for problems with delay derivative order > 1, try increasing delaypts", evalf[8](_val) elif _dat[4][9] = 10 then error "cannot evaluate the solution further right of %1, interrupt requested", evalf[8](_val) elif 100 < _dat[4][9] then if _dat[4][9]-100 = _nv+1 then error "constraint projection failure on event at t=%1", evalf[8](_val) elif _dat[4][9]-100 = _nv+2 then error "index-1 and derivative evaluation failure on event at t=%1", evalf[8](_val) elif _dat[4][9]-100 = _nv+3 then error "maximum number of event iterations reached (%1) at t=%2", round(_dat[3][1][_nv+1, 3]), evalf[8](_val) else if _Env_dsolve_nowarnstop <> true then `dsolve/numeric/warning`(StringTools:-FormatMessage("cannot evaluate the solution further right of %1, event #%2 triggered a halt", evalf[8](_val), round(_dat[3][1][_dat[4][9]-100, 1]))) end if; Rounding := 'nearest'; _xout := _val end if else error "cannot evaluate the solution further right of %1", evalf[8](_val) end if elif _i = 2 and _xout < _val then Rounding := infinity; if _dat[4][9] = 1 then error "cannot evaluate the solution further left of %1, probably a singularity", evalf[8](_val) elif _dat[4][9] = 2 then error "cannot evaluate the solution further left of %1, maxfun limit exceeded (see ?dsolve,maxfun for details)", evalf[8](_val) elif _dat[4][9] = 3 then if _dat[4][25] = 3 then error "cannot evaluate the solution past the initial point, problem may be initially singular or improperly set up" else error "cannot evaluate the solution past the initial point, problem may be complex, initially singular or improperly set up" end if elif _dat[4][9] = 4 then error "cannot evaluate the solution further left of %1, accuracy goal cannot be achieved with specified 'minstep'", evalf[8](_val) elif _dat[4][9] = 5 then error "cannot evaluate the solution further left of %1, too many step failures, tolerances may be too loose for problem", evalf[8](_val) elif _dat[4][9] = 6 then error "cannot evaluate the solution further left of %1, cannot downgrade delay storage for problems with delay derivative order > 1, try increasing delaypts", evalf[8](_val) elif _dat[4][9] = 10 then error "cannot evaluate the solution further right of %1, interrupt requested", evalf[8](_val) elif 100 < _dat[4][9] then if _dat[4][9]-100 = _nv+1 then error "constraint projection failure on event at t=%1", evalf[8](_val) elif _dat[4][9]-100 = _nv+2 then error "index-1 and derivative evaluation failure on event at t=%1", evalf[8](_val) elif _dat[4][9]-100 = _nv+3 then error "maximum number of event iterations reached (%1) at t=%2", round(_dat[3][1][_nv+1, 3]), evalf[8](_val) else if _Env_dsolve_nowarnstop <> true then `dsolve/numeric/warning`(StringTools:-FormatMessage("cannot evaluate the solution further left of %1, event #%2 triggered a halt", evalf[8](_val), round(_dat[3][1][_dat[4][9]-100, 1]))) end if; Rounding := 'nearest'; _xout := _val end if else error "cannot evaluate the solution further left of %1", evalf[8](_val) end if end if; if _EnvInFsolve = true then _dig := _dat[4][26]; _dat[4][26] := _EnvDSNumericSaveDigits; _Env_dsolve_SC_native := true; if _dat[4][25] = 1 then _i := 1; _dat[4][25] := 2 else _i := _dat[4][25] end if; _val := `dsolve/numeric/SC/IVPval`(_dat, _xout, _src); _dat[4][25] := _i; _dat[4][26] := _dig; [_xout, seq(_val[_vmap[_i]], _i = 1 .. _n-_ne)] else Digits := _dat[4][26]; _val := `dsolve/numeric/SC/IVPval`(eval(_dat, 2), _xout, _src); [_xout, seq(_val[_vmap[_i]], _i = 1 .. _n-_ne)] end if end proc, (2) = Array(0..0, {}), (3) = [t, x(t)], (4) = []}); _vars := _dat[3]; _pars := map(rhs, _dat[4]); _n := nops(_vars)-1; _solnproc := _dat[1]; if not type(_xout, 'numeric') then if member(x_rkf45, ["start", 'start', "method", 'method', "left", 'left', "right", 'right', "leftdata", "rightdata", "enginedata", "eventstop", 'eventstop', "eventclear", 'eventclear', "eventstatus", 'eventstatus', "eventcount", 'eventcount', "laxtol", 'laxtol', "numfun", 'numfun', NULL]) then _res := _solnproc(convert(x_rkf45, 'string')); if 1 < nops([_res]) then return _res elif type(_res, 'array') then return eval(_res, 1) elif _res <> "procname" then return _res end if elif member(x_rkf45, ["last", 'last', "initial", 'initial', "parameters", 'parameters', "initial_and_parameters", 'initial_and_parameters', NULL]) then _xout := convert(x_rkf45, 'string'); _res := _solnproc(_xout); if _xout = "parameters" then return [seq(_pars[_i] = _res[_i], _i = 1 .. nops(_pars))] elif _xout = "initial_and_parameters" then return [seq(_vars[_i+1] = [_res][1][_i+1], _i = 0 .. _n), seq(_pars[_i] = [_res][2][_i], _i = 1 .. nops(_pars))] else return [seq(_vars[_i+1] = _res[_i+1], _i = 0 .. _n)] end if elif type(_xout, `=`) and member(lhs(_xout), ["initial", 'initial', "parameters", 'parameters', "initial_and_parameters", 'initial_and_parameters', NULL]) then _xout := convert(lhs(x_rkf45), 'string') = rhs(x_rkf45); if type(rhs(_xout), 'list') then _res := _solnproc(_xout) else error "initial and/or parameter values must be specified in a list" end if; if lhs(_xout) = "initial" then return [seq(_vars[_i+1] = _res[_i+1], _i = 0 .. _n)] elif lhs(_xout) = "parameters" then return [seq(_pars[_i] = _res[_i], _i = 1 .. nops(_pars))] else return [seq(_vars[_i+1] = [_res][1][_i+1], _i = 0 .. _n), seq(_pars[_i] = [_res][2][_i], _i = 1 .. nops(_pars))] end if elif type(_xout, `=`) and member(lhs(_xout), ["eventdisable", 'eventdisable', "eventenable", 'eventenable', "eventfired", 'eventfired', "direction", 'direction', NULL]) then return _solnproc(convert(lhs(x_rkf45), 'string') = rhs(x_rkf45)) elif _xout = "solnprocedure" then return eval(_solnproc) elif _xout = "sysvars" then return _vars end if; if procname <> unknown then return ('procname')(x_rkf45) else _ndsol; _ndsol := pointto(_dat[2][0]); return ('_ndsol')(x_rkf45) end if end if; try _res := _solnproc(_xout); [seq(_vars[_i+1] = _res[_i+1], _i = 0 .. _n)] catch: error  end try end proc

FileTools:-Text:-ReadFile(str);

"proc(x_rkf45)
local _res, _dat, _vars, _solnproc, _xout, _ndsol, _pars, _n, _i;
option `Copyright (c) 2000 by Waterloo Maple Inc. All rights reserved.`;
  if 1 < nargs then error "invalid input: too many arguments" end if;
  _EnvDSNumericSaveDigits := Digits;
  Digits := 15;
  if _EnvInFsolve = true then
    _xout := evalf[_EnvDSNumericSaveDigits](x_rkf45)
  else
    _xout := evalf(x_rkf45)
  end if;
  _dat := Array(1..4, {(1) = proc (_xin) local _xout, _dtbl, _dat, _vmap, 
   _x0, _y0, _val, _dig, _n, _ne, _nd, _nv, _pars, _ini, _par, _i, _j, _k, 
   _src; option `Copyright (c) 2002 by Waterloo Maple Inc. All rights 
   reserved.`; table( [( "complex" ) = false ] ) _xout := _xin; _pars := []; 
   _dtbl := array( 1 .. 4, [( 1 ) = (array( 1 .. 24, [( 1 ) = (datatype = 
   float[8], order = C_order, storage = rectangular), ( 2 ) = (datatype = 
   float[8], order = C_order, storage = rectangular), ( 3 ) = ([0, 0, 0, 
   Array(1..0, 1..2, {}, datatype = float[8], order = C_order)]), ( 4 ) = 
   (Array(1..54, {(1) = 1, (2) = 1, (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) = 21, (19) = 30000, (20) = 0, (21) = 
   0, (22) = 1, (23) = 4, (24) = 0, (25) = 1, (26) = 15, (27) = 1, (28) = 0, 
   (29) = 1, (30) = 3, (31) = 3, (32) = 0, (33) = 1, (34) = 0, (35) = 0, 
   (36) = 0, (37) = 0, (38) = 0, (39) = 0, (40) = 0, (41) = 0, (42) = 0, 
   (43) = 1, (44) = 0, (45) = 0, (46) = 0, (47) = 0, (48) = 0, (49) = 0, 
   (50) = 50, (51) = 1, (52) = 0, (53) = 0, (54) = 0}, datatype = 
   integer[8])), ( 5 ) = (Array(1..28, {(1) = .0, (2) = 0.10e-5, (3) = .0, 
   (4) = 0.500001e-14, (5) = .0, (6) = .16695651424556227, (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..1, {(1) = -1.0}, 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.318908691406\
  25e-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.2212245092262\
  748, (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.999903528199\
  06, (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.992771707568\
  7275, (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..1, {(1) = .1}, datatype = float[8], order = C_order), 
   Array(1..1, {(1) = .0}, datatype = float[8], order = C_order), 
   Array(1..1, {(1) = .0}, datatype = float[8], order = C_order), 
   Array(1..1, {(1) = .0}, datatype = float[8], order = C_order), 
   Array(1..1, {(1) = .0}, datatype = float[8], order = C_order), 
   Array(1..1, 1..1, {(1, 1) = .0}, datatype = float[8], order = C_order), 
   Array(1..1, 1..1, {(1, 1) = .0}, datatype = float[8], order = C_order), 
   Array(1..1, 1..6, {(1, 1) = .0, (1, 2) = .0, (1, 3) = .0, (1, 4) = .0, 
   (1, 5) = .0, (1, 6) = .0}, datatype = float[8], order = C_order), 
   Array(1..1, {(1) = 0}, datatype = integer[8]), Array(1..1, {(1) = .0}, 
   datatype = float[8], order = C_order), Array(1..1, {(1) = .0}, datatype = 
   float[8], order = C_order), Array(1..1, {(1) = .0}, datatype = float[8], 
   order = C_order), Array(1..1, {(1) = .0}, datatype = float[8], order = 
   C_order), Array(1..1, {(1) = .0}, datatype = float[8], order = 
   C_order)]), ( 8 ) = ([Array(1..1, {(1) = -1.0}, datatype = float[8], 
   order = C_order), Array(1..1, {(1) = .0}, datatype = float[8], order = 
   C_order), Array(1..1, {(1) = .0}, datatype = float[8], order = C_order), 
   0, 0]), ( 11 ) = (Array(1..6, 0..1, {(1, 1) = .0, (2, 0) = .0, (2, 1) = 
   .0, (3, 0) = .0, (3, 1) = .0, (4, 0) = .0, (4, 1) = .0, (5, 0) = .0, (5, 
   1) = .0, (6, 0) = .0, (6, 1) = .0}, datatype = float[8], order = 
   C_order)), ( 10 ) = ([proc (N, X, Y, YP) option `[Y[1] = x(t)]`; YP[1] := 
   sin(X); 0 end proc, -1, 0, 0, 0, 0, 0, 0]), ( 13 ) = (), ( 12 ) = (), ( 
   15 ) = ("rkf45"), ( 14 ) = ([0, 0]), ( 18 ) = ([]), ( 19 ) = (0), ( 16 ) 
   = ([0, 0, 0, []]), ( 17 ) = ([proc (N, X, Y, YP) option `[Y[1] = x(t)]`; 
   YP[1] := sin(X); 0 end proc, -1, 0, 0, 0, 0, 0, 0]), ( 22 ) = (0), ( 23 ) 
   = (0), ( 20 ) = ([]), ( 21 ) = (0), ( 24 ) = (0)  ] ))  ] ); _y0 := 
   Array(0..1, {(1) = 0.}); _vmap := array( 1 .. 1, [( 1 ) = (1)  ] ); _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 `or`(_Env_smart_\
  dsolve_numeric = true, _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 `or`(`and`(`<>`(_dtbl[4], 2), `<>`(_dtbl[4], 3)), `+`(_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 `and`(`and`(type(_xin, `=`), 
   type(rhs(_xin), 'list')), member(lhs(_xin), {"initial", "parameters", 
   "initial_and_parameters"})) then _ini, _par := [], []; if lhs(_xin) = 
   "initial" then _ini := rhs(_xin) elif lhs(_xin) = "parameters" then _par 
   := rhs(_xin) elif `<>`(select(type, rhs(_xin), `=`), []) then _par, _ini 
   := selectremove(type, rhs(_xin), `=`) elif `<`(nops(rhs(_xin)), 
   `+`(nops(_pars), 1)) then error "insufficient data for specification of 
   initial and parameters" else _par := rhs(_xin)[`+`(`-`(nops(_pars))) .. 
   -1]; _ini := rhs(_xin)[1 .. `+`(`-`(nops(_pars)), `-`(1))] end if; _xout 
   := lhs(_xout); if `<>`(_par, []) then `dsolve/numeric/process_parameters`(\
  _n, _pars, _par, _y0) end if; if `<>`(_ini, []) then `dsolve/numeric/proces\
  s_initial`(`+`(_n, `-`(_ne)), _ini, _y0, _pars, _vmap) end if; 
   `dsolve/numeric/SC/reinitialize`(_dtbl, _y0, _n, procname, _pars); if 
   `and`(`and`(_Env_smart_dsolve_numeric = true, type(_y0[0], 'numeric')), 
   `<>`(_dtbl[1][4][10], 1)) then procname("right") := _y0[0]; procname("left\
  ") := _y0[0] end if; if _xout = "initial" then return [_y0[0], 
   seq(_y0[_vmap[_i]], _i = 1 .. `+`(_n, `-`(_ne)))] elif _xout = 
   "parameters" then return [seq(_y0[`+`(_n, _i)], _i = 1 .. nops(_pars))] 
   else return [_y0[0], seq(_y0[_vmap[_i]], _i = 1 .. `+`(_n, `-`(_ne)))], 
   [seq(_y0[`+`(_n, _i)], _i = 1 .. nops(_pars))] end if elif _xin = 
   "eventstop" then if _nv = 0 then error "this solution has no events" end 
   if; _i := _dtbl[4]; if `and`(`<>`(_i, 2), `<>`(_i, 3)) then return 0 end 
   if; if `and`(`and`(`and`(_dtbl[_i][4][10] = 1, assigned(_dtbl[`+`(5, 
   `-`(_i))])), `<`(_dtbl[_i][4][9], 100)), `<=`(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 `and`(`<>`(_i, 2), `<>`(_i, 3)) then 
   error "no events to clear" end if; if `and`(`and`(`and`(_dtbl[_i][4][10] 
   = 1, assigned(_dtbl[`+`(5, `-`(_i))])), `<`(_dtbl[_i][4][9], 100)), 
   `<`(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 `and`(type(_xin, `=`), 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(_dt\
  bl[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(`and`(a\
  ssigned(_dtbl[2]), member(_dtbl[2][4][17], _i))), evalb(`and`(assigned(_dtb\
  l[3]), 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 `and`(`<=`(_k, _nv), 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 `and`(_dtbl[2][3][1][_k, 2] = 0, 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 `and`(_dtbl[2][3][1][_k, 2] 
   = 0, 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 `and`(`<=`(_k, _nv), 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 `and`(_dtbl[3][3][1][_\
  k, 2] = 0, 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 `and`(_dtbl[3][3][1][_k, 2] = 0, irem(iquo(round(_dtb\
  l[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 `and`(type(_xin, `=`), 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 `and`(`<>`(_dtbl[4], 2), `<>`(_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 `and`(type(_k, 'integer' = 
   'anything'), 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 `or`(`<`(_src, 1), `<`(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 `and`(`<>`(_dtbl[1][3][1][_src, 2], 0.), 
   `<>`(`+`(_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 _EnvEventRetrigge\
  rWarned = false then WARNING(`'eventfired' has no effect on events that 
   retrigger`) end if; _EnvEventRetriggerWarned := true end if; if 
   `and`(_dtbl[_i][3][1][_src, 2] = 0, irem(iquo(round(_dtbl[_i][3][1][_src, 
   4]), 32), 2) = 1) then _val := _val, undefined elif `or`(`or`(type(_dtbl[_\
  i][3][4][_src, `+`(_i, `-`(1))], 'undefined'), `and`(_i = 2, `<`(_dtbl[2][3\
  ][1][_src, 8], _dtbl[2][3][4][_src, 1]))), `and`(_i = 3, `<`(_dtbl[3][3][4]\
  [_src, 2], _dtbl[3][3][1][_src, 8]))) then _val := _val, _dtbl[_i][3][1][_s\
  rc, 8] else _val := _val, _dtbl[_i][3][4][_src, `+`(_i, `-`(1))] end if; 
   if type(_k, `=`) then if `and`(_dtbl[_i][3][1][_src, 2] = 0, 
   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 
   `and`(type(_xin, `=`), lhs(_xin) = "direction") then if not member(rhs(_xi\
  n), {-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 `and`(`<=`(_j, _nv), 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 `and`(_dtbl[_i][3][1][_j, 2] = 0, 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 `and`(_dtbl[_i][3][1][_\
  j, 2] = 0, 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 `or`(_dtbl[1][3][1] = 0, `and`(`<>`(_dtbl\
  [4], 2), `<>`(_dtbl[4], 3))) then return 0 else return round(_dtbl[_dtbl[4]\
  ][3][1][`+`(_nv, 1), 12]) end if else return "procname" end if end if; if 
   _xout = _x0 then return [_x0, seq(evalf(_dtbl[1][6][_vmap[_i]]), _i = 1 
   .. `+`(_n, `-`(_ne)))] end if; _i := `if`(`<=`(_x0, _xout), 3, 2); if 
   `and`(`and`(_xin = "last", `<`(0, _dtbl[_i][4][9])), `<`(_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 `and`(`<=`(_j, _nv), 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 `and`(_dtbl[_i][3][1][_j, 2] = 0, 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 `and`(_dtbl[_i][3][1][_\
  j, 2] = 0, 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`(_d\
  tbl, _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 `and`(_src = 0, `<`(100, _dat[4][9])) then _val 
   := _dat[3][1][`+`(_nv, 1), 8] else _val := _dat[11][_dat[4][20], 0] end 
   if; if `or`(`<>`(_src, 0), `<=`(_dat[4][9], 0)) then _dtbl[1][5][1] := 
   _xout else _dtbl[1][5][1] := _val end if; if `and`(_i = 3, `<`(_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 `and`(_i = 2, `<`(_xout, 
   _val)) then Rounding := infinity; if _dat[4][9] = 1 then error "cannot 
   evaluate the solution further left of %1, probably a singularity", 
   evalf[8](_val) elif _dat[4][9] = 2 then error "cannot evaluate the 
   solution further left of %1, maxfun limit exceeded (see ?dsolve,maxfun 
   for details)", evalf[8](_val) elif _dat[4][9] = 3 then if _dat[4][25] = 3 
   then error "cannot evaluate the solution past the initial point, problem 
   may be initially singular or improperly set up" else error "cannot 
   evaluate the solution past the initial point, problem may be complex, 
   initially singular or improperly set up" end if elif _dat[4][9] = 4 then 
   error "cannot evaluate the solution further left of %1, accuracy goal 
   cannot be achieved with specified 'minstep'", evalf[8](_val) elif 
   _dat[4][9] = 5 then error "cannot evaluate the solution further left of 
   %1, too many step failures, tolerances may be too loose for problem", 
   evalf[8](_val) elif _dat[4][9] = 6 then error "cannot evaluate the 
   solution further left of %1, cannot downgrade delay storage for problems 
   with delay derivative order > 1, try increasing delaypts", evalf[8](_val) 
   elif _dat[4][9] = 10 then error "cannot evaluate the solution further 
   right of %1, interrupt requested", evalf[8](_val) elif `<`(100, 
   _dat[4][9]) then if `+`(_dat[4][9], `-`(100)) = `+`(_nv, 1) then error 
   "constraint projection failure on event at t=%1", evalf[8](_val) elif 
   `+`(_dat[4][9], `-`(100)) = `+`(_nv, 2) then error "index-1 and 
   derivative evaluation failure on event at t=%1", evalf[8](_val) elif 
   `+`(_dat[4][9], `-`(100)) = `+`(_nv, 3) then error "maximum number of 
   event iterations reached (%1) at t=%2", round(_dat[3][1][`+`(_nv, 1), 
   3]), evalf[8](_val) else if `<>`(_Env_dsolve_nowarnstop, true) then 
   `dsolve/numeric/warning`(StringTools:-FormatMessage("cannot evaluate the 
   solution further left of %1, event #%2 triggered a halt", evalf[8](_val), 
   round(_dat[3][1][`+`(_dat[4][9], `-`(100)), 1]))) end if; Rounding := 
   'nearest'; _xout := _val end if else error "cannot evaluate the solution 
   further left of %1", evalf[8](_val) end if end if; if _EnvInFsolve = true 
   then _dig := _dat[4][26]; _dat[4][26] := _EnvDSNumericSaveDigits; 
   _Env_dsolve_SC_native := true; if _dat[4][25] = 1 then _i := 1; 
   _dat[4][25] := 2 else _i := _dat[4][25] end if; _val := `dsolve/numeric/SC\
  /IVPval`(_dat, _xout, _src); _dat[4][25] := _i; _dat[4][26] := _dig; 
   [_xout, seq(_val[_vmap[_i]], _i = 1 .. `+`(_n, `-`(_ne)))] else Digits := 
   _dat[4][26]; _val := `dsolve/numeric/SC/IVPval`(eval(_dat, 2), _xout, 
   _src); [_xout, seq(_val[_vmap[_i]], _i = 1 .. `+`(_n, `-`(_ne)))] end if 
   end proc, (2) = Array(0..0, {}), (3) = [t, x(t)], (4) = []});
  _vars := _dat[3];
  _pars := map(rhs, _dat[4]);
  _n := nops(_vars) - 1;
  _solnproc := _dat[1];
  if not type(_xout, 'numeric') then
    if member(x_rkf45, ["start", 'start', "method", 'method', "left", 'left', 
      "right", 'right', "leftdata", "rightdata", "enginedata", "eventstop", 
      'eventstop', "eventclear", 'eventclear', "eventstatus", 'eventstatus', 
      "eventcount", 'eventcount', "laxtol", 'laxtol', "numfun", 'numfun', NULL
      ]) then
      _res := _solnproc(convert(x_rkf45, 'string'));
      if 1 < nops([_res]) then
        return _res
      elif type(_res, 'array') then
        return eval(_res, 1)
      elif _res <> "procname" then
        return _res
      end if;
    elif member(x_rkf45, ["last", 'last', "initial", 'initial', "parameters", 
      'parameters', "initial_and_parameters", 'initial_and_parameters', NULL]) then
      _xout := convert(x_rkf45, 'string');
      _res := _solnproc(_xout);
      if _xout = "parameters" then
        return [seq(_pars[_i] = _res[_i], _i = 1 .. nops(_pars))]
      elif _xout = "initial_and_parameters" then
        return [seq(_vars[_i + 1] = [_res][1][_i + 1], _i = 0 .. _n), 
        seq(_pars[_i] = [_res][2][_i], _i = 1 .. nops(_pars))]
      else
        return [seq(_vars[_i + 1] = _res[_i + 1], _i = 0 .. _n)]
      end if;
    elif type(_xout, `=`) and member(lhs(_xout), ["initial", 'initial', 
      "parameters", 'parameters', "initial_and_parameters", 
      'initial_and_parameters', NULL]) then
      _xout := convert(lhs(x_rkf45), 'string') = rhs(x_rkf45);
      if type(rhs(_xout), 'list') then
        _res := _solnproc(_xout)
      else
        error "initial and/or parameter values must be specified in a list"
      end if;
      if lhs(_xout) = "initial" then
        return [seq(_vars[_i + 1] = _res[_i + 1], _i = 0 .. _n)]
      elif lhs(_xout) = "parameters" then
        return [seq(_pars[_i] = _res[_i], _i = 1 .. nops(_pars))]
      else
        return [seq(_vars[_i + 1] = [_res][1][_i + 1], _i = 0 .. _n), 
        seq(_pars[_i] = [_res][2][_i], _i = 1 .. nops(_pars))]
      end if;
    elif type(_xout, `=`) and member(lhs(_xout), ["eventdisable", 
      'eventdisable', "eventenable", 'eventenable', "eventfired", 
      'eventfired', "direction", 'direction', NULL]) then
      return _solnproc(convert(lhs(x_rkf45), 'string') = rhs(x_rkf45))
    elif _xout = "solnprocedure" then
      return eval(_solnproc)
    elif _xout = "sysvars" then
      return _vars
    end if;
    if procname <> unknown then
      return 'procname'(x_rkf45)
    else
      _ndsol;
      _ndsol := pointto(_dat[2][0]);
      return '_ndsol'(x_rkf45);
    end if;
  end if;
  try
    _res := _solnproc(_xout);
    [seq(_vars[_i + 1] = _res[_i + 1], _i = 0 .. _n)];
  catch:
    error 
  end try;
end proc;


"

 


 

Download writetoexample.mw

@DSkoog I think that the explanation you've given about why it does not work is not right. For example, try using `printf` after calling `writeto` to redirect. The `printf` command does not produce "assignable returned output", but its output can be redirected.

No, the problem is rather that `showstat` makes use of an internal routine that calls `fprintf` with an (essentially) preset name. But this can be undone. I'll put one way to do it in an Answer.

Of course, it's still possible that the OP doesn't especially want or need line numbers, and simply using `print` (after setting verboseproc) instead of `showstat` might well be adequate for him.

I changed the inlined URLs to the original authors' webpages, since there is no evidence that the poster has those authors' permission to copy any of the materials up to this site.

Applying the evalf command to such a RootOf will cause Maple to use fsolve.  (This doesn't you all control over the various options to the fsolve command, though.)

Another way is to fix the hue layer of the color Array in the plot structure, after the fact.
 

restart;

a:=0: b:=1:

with(plots):

P:=densityplot(z,dummy=0..1,z=a..b,grid=[2,10],size=[90,100],
               colorscheme=["zcoloring",[z->(a-z)/(b-a),
                                         z->1,z->1],
                            colorspace = "HSV"],
               style=surface,axes=frame,labels=["",""],
               axis[1]=[tickmarks=[]],title="Test",
               titlefont = ["Calibri","Bold",16],
               axesfont=["Arial",14],size=[100,500],
               restricttoranges):

newdat := copy(op([1,4,2],P)):

newdat[..,..,1]:=240/360*newdat[..,..,1]:

newP := subsop([1,4,2]=newdat,P);

 


 

Download coolhot2.mw

@Rouben Rostamian  Since the OP appears to have Maple 18, then the commands would be,

map(fnormal, M.v);

simplify(map(fnormal, M.v), zero);

map(fnormal, v);

simplify(map(fnormal, v), zero);

I'm not sure whether this might confuse the OP, but the smallest eigenvalue (in absolute value) of the given Matrix is not zero. It is small in absolute value, but not zero. So it's not clear whether he would prefer to find the eigenvector associated with that, or the singular vector (which is what NullSpace computes here). I note also that for higher working precision (Digits>12 ?) the numerical rank will compute as 4 and NullSpace will return empty.

@John Fredsted Fooling around and learning by experimentation is always good. Nobody was born knowing maple.

My attachment won't serve as oracle for other different examples. Perhaps the key thing to remember is that the commands zip, map (and map [n] variants) and elementwise operators are different. There is often overlap in their functionality, but no simple correspondence. 

@John Fredsted It seems to me that quite a bit of your difficulty in this example is due to misconceptions about what elementwise operators are supposed to do, and what the syntax means.

Your original map call was,

    map(x -> A . x,LM)

and one valid elementwise operator equivalent is,

    (x -> A . x)~(LM)

since the original operator in your question is x->A.x and not the `.` operator.

It may be that, conceptually, you were only considering that LM is a "container", and overlooking that it really does matter that A is also a "container". But that matters, in the attempted A .~ LM you tried. Bear in mind that your goal was a computation which does not act elementwise on A, so plain A would not be an argument to the successful elementwise operator call. Your attempt A .~ LM is the infix variant of the prefix call `.`~(A,LM) and so it's quite logical that would attempt to use both A and LM in an elementwise manner.

As for the map2 example I gave, it might help a little to think of map as map[1], and map2 as map[2]. You wanted to treat only LM (and not A) in an elementwise manner. The map[1] and map[2] commands allow you to have just one of the arguments be treated elementwise. And that extends to map[n] more generally,

    `*`~( [2,3], [a,b], [x,y] );

                               [2 a x, 3 b y]

    map[1](`*`, [2,3], [a,b], [x,y] );

                     [2 [a, b] [x, y], 3 [a, b] [x, y]]

    map[2](`*`, [2,3], [a,b], [x,y] );

                     [[2, 3] a [x, y], [2, 3] b [x, y]]

    map[3](`*`, [2,3], [a,b], [x,y] );

                     [[2, 3] [a, b] x, [2, 3] [a, b] y]

I'll attach a sheet. I hope it helps clarify, rather than not.

elementwise.mw

@tomleslie On the help page for topic operators,elementwise in Maple 2016 I see, "When mixed container types are present the return type will be determined according to the following precedence: object, rtable, list, set."

And that explains why both [1,2,3]*~{4,5,6} and {4,5,6}*~[1,2,3] return a list. It also (correctly) indicates that elementwise multiplication of a Vector and a set returns a Vector.

First 287 288 289 290 291 292 293 Last Page 289 of 592