Events, in combination with discrete variables, are a mechanism that can be used in dsolve/numeric to handle stopping criteria, reset conditions, zero order holds, and most other events that occur in ODE and DAE system simulation. Note that events are restricted to operations that act on variable values, and cannot be used to change the essential form of the ODE system or exchange equations.

These are available in dsolve,numeric for the default non-stiff (rkf45, rkf45_dae, ck45, ck45_dae) and stiff (rosenbrock and rosenbrock_dae) solvers. Note that events replace the stop condition mechanism previously available for these solvers.

The sections below describe the behavior and features of the event implementation in dsolve[numeric].

Events are specified by the events=[event1,event2,...] option in the call to dsolve[numeric], where each event is a list of [trigger,action] pairs, where the trigger describes the trigger of the event, and the action describes the action to perform when the event is triggered.

The following are the options for the form of the trigger part of the description:

y(t): simple root finding trigger

f(t,y(t)): root finding trigger

[f(t,y(t)),c(t,y(t))<0]: root finding with a conditional trigger

[f(t,y(t)), And(c1(t,y(t)) > 0, c2(t,y(t)) > 0)]: root finding with a compound conditional trigger

[0, c(t,y(t)) < 0]: discrete event with a conditional trigger

[0, And(c1(t,y(t)) > 0, c2(t,y(t)) > 0)]: discrete event with a compound conditional trigger

f(t,y(t)) = lo..hi: range trigger, that fires once the expression goes outside the range

t=start: timed trigger that fires when t=start

t=[start, $\mathrm{\delta }$]: timed trigger that fires when t=..., start-$\mathrm{\delta }$, start, start+$\mathrm{\delta }$, start+2*$\mathrm{\delta }$, ...

Any events where the trigger depends on continuous variables fire when the trigger becomes satisfied during the integration process (on the so-called leading edge). This means that if a continuous condition is satisfied at the initial point, it will not fire (not by default, though this can be changed via fireinitial described below), but will fire if the condition is satisfied at any time past the initial point.

This is not the case for timed triggers, which do fire at the initial point, so, for example, the trigger t-t0 will not fire at the initial point, while the trigger t=t0 will.

Discrete events (where either the rootfind part is zero, or the event contains only discrete variables) can fire only during an event iteration that has been initiated by a continuous trigger (see the following sections), or if it has been specified that events iteration should occur each time the integrator resumes (the event_initial option).

The following are the options for the form of the action part of the description:

none: do nothing, but this forces a point to be included at the trigger point for an adaptive plot mode (for example, the range option of adaptive odeplot output), and induces an event iteration which can trigger discrete events.

halt: halt integration with a stop message when triggered

tobegin: (iteration) jump to the start of the event iteration, incrementing the iteration counter and updating pre.

toend: (iteration) jump to the end of the event iteration, performing a termination check before iterating.

breakiteration: (iteration) halt event iteration, resuming continuous integration.

delayhalt: (iteration) schedule a halt of the integration on completion of event iteration.

nomark: (iteration) do not tag event as having executed in an event iteration (only useful if event_iterate=trigger and described in more detail in the discussion of event iteration).

u(t) = -u(t): execute statement when triggered

[u(t) = -u(t), y(t) = -y(t)]: execute statements when triggered

[If(y(t) < 0, y(t) = -y(t), halt)]: execute statements conditionally when triggered

The conditional part of a trigger can use the inert (capitalized) form of the Not, And, Or, Xor, and Implies logical operations. The action can use these, the If conditional, the pre operator (see below), halt, tobegin, toend, breakiteration, delayhalt, nomark and none; within the action, equations are interpreted as assignments. Note that delayhalt and nomark must be the last statement in any embedded code, as when executed they force an immediate return to the controlling routine (which is obviously also true for halt, tobegin, toend, and breakiteration). Note that only none cannot be embedded in code, and only nochange must be embedded in code. Note that the action can also use temporary variables in the associated computation, so long as the names are distinct from the dependent and independent variable names for the problem. It can also use derivative values up to one order lower than the order of the derivative in the problem (that is, for a second order ODE in y(t), up to diff(y(t), t) can be used). As an example, a valid action for a second order ODE in y(t) could be [t1 = sqrt(y(x)^2 + diff(y(t),t)^2), x(t) = x(t)*t1, y(t) = y(t)*t1].

When a root-finding trigger fires, the value of the expression in the root finder is only an approximate zero due to numerical round-off. What this means is that f(t,y(t)) may evaluate, with respect to the current solution, to a value slightly above or slightly below zero. This may not be desirable for some computations where it is critical that the value of the trigger expression is either exactly zero, as approaching from the left, or as just past the trigger point. This is primarily desired to be able to apply events for an ODE system that has been separated into disjoint cases dependent on the values of particular triggers (in which case you always want to use a form that provides the values just past the trigger point).

A mechanism has been provided for this, but only for simple root finding triggers (i.e. the expression for the root find must be a dependent variable of the problem). In the event that an expression of the dependent and independent variables is required, an index-1 relation of the form newvar(t) = f(t,y(t)) can be added, and the trigger can then be based on newvar(t).

The side function, in place of the root finding dependent variable y(t) can be used in the forms side(y(t),pre) or side(y(t),post) to obtain the value just prior, at, or just after the root.

For example, for a root finding trigger side(y(t),post) where y(t) is initially negative, at the time at which the event is triggered the value of y(t) will be a small positive quantity.

Discrete variables are a way of maintaining discrete data as part of a numerical ODE/DAE computation in dsolve[numeric], and can be as simple as a boolean flag that indicates that a switch is either in the on or off position, an integer value that counts how many times a particular criterion is met, or a sampling of continuous data, as would be required for a zero-order hold.

The values of discrete variables can only ever be changed in events, since they are not continuous (and hence will not have an associated differential equation). As a result, derivatives of discrete variables cannot appear in the input to dsolve[numeric].

Discrete variables depend upon the independent variable (since their value can change within event handling), so must be declared as functions of the independent variable. These can be specified to dsolve[numeric] by using the option discrete_variables=[v1,v2,...], where for any i, vi can be a dependent variable (for example, u(t)) or a dependent variable with a declared type (u(t)::integer). Only the types boolean, integer and float are currently available. Note: it is strongly recommended that any discrete variables that are not of type float (the default) be explicitly declared.

Discrete variables must also be given an initial value (even if it is not needed) in the input system.

As an example of the use of a discrete variable for a zero order hold on the continuous dependent variable y(t) sampled at $\frac{1}{2}$ second intervals, specify discrete_variables=[u(t)::float], events=[[t = [0,1/2], u(t) = y(t)]] and an initial condition on u(t).

Note that discrete variables that are declared as boolean can be used directly in the conditionals of any event construct (for example, if u(t) was declared as boolean, then a conditional can be of the form If(u(t),...)).

Another mechanism for handling events that is most useful in cases where a trigger expression asymptotically approaches (i.e. touches) zero is the ability to associate a discrete (boolean) variable to a trigger.

This is done using the requires option, and is only available for simple triggers (i.e. basic root finding triggers, supporting increasing, but with no additional options, and no trigger condition). Also note that increasing triggers must occur in pairs, where one trigger depends on the requires variable, and the other depends on the Not of the requires variable, and the sum of the arguments to the corresponding increasing (or the difference of the increasing and decreasing arguments) must be a constant.

The value of the requires option specifies the variable associated to the trigger, and the value that the variable should have when the trigger fires in an increasing direction (which is the opposite to the value the variable will have after the firing of the trigger). The requires variable will be automatically updated to the opposite value as specified for the trigger firing prior to execution of the corresponding action, and the action cannot update the value of the requires variable (though it can use it).

As a simple example, consider a problem with continuous dependent variable $y\left(t\right)$, where the discrete variable $\mathrm{py}\left(t\right)$ is true (1) if $0<y\left(t\right)$ and false (0) if $y\left(t\right)<0$ (and note that the value when $y\left(t\right)=0$ can be either). The event specification that can be used to tie $\mathrm{py}\left(t\right)$ into $y\left(t\right)$ is simply $\left[y\left(t\right)\&comma;\mathrm{none}\&comma;\mathrm{requires}=\left(\neg \mathrm{py}\left(t\right)\right)\right]$, and the problem must start with an initial condition that is consistent (i.e. if $0<y\left(t\right)$ then $\mathrm{py}\left(t\right)=1$). This trigger will then fire when $y\left(t\right)$ passes through zero in the positive direction with $\mathrm{py}\left(t\right)=0$, and will also fire when $y\left(t\right)$ passes through zero in the negative direction with $\mathrm{py}\left(t\right)=1$.

The same behavior can be accomplished with a somewhat less efficient pair of events: $\left[\mathrm{increasing}\left(y\left(t\right)\right)\&comma;\mathrm{requires}=\left(\neg \mathrm{py}\left(t\right)\right)\right]$ and $\left[\mathrm{increasing}\left(-y\left(t\right)\right)\&comma;\mathrm{requires}=\mathrm{py}\left(t\right)\right]$. The first event sets $\mathrm{py}\left(t\right)=1$ when $y\left(t\right)$ passes through zero in the positive direction with $\mathrm{py}\left(t\right)=0$, and the second event sets $\mathrm{py}\left(t\right)=0$ when $y\left(t\right)$ passes through zero in the negative direction with $\mathrm{py}\left(t\right)=1$.

The latter mechanism is most useful when providing a band to avoid spurious events triggering as a result of round off (so one would change the trigger expressions to $\mathrm{increasing}\left(y\left(t\right)-\mathrm{\epsilon }\right)$ and $\mathrm{increasing}\left(-y\left(t\right)-\mathrm{\epsilon }\right)$ respectively to provide an $\mathrm{\epsilon }$ width band about the trigger expression).

Events that are based on continuous criteria (such as root finding or events that fire for a fixed time or sequence of times) can fire only once at any given time. This is not the case for events that trigger on discrete variables.

Events that are based on discrete variables can be utilized to code an algorithm (of sorts) into the event process, but note that for the events to be initiated, the condition must be satisfied at a point after the trigger of a continuous event (for example, a sign change or a certain time being reached).

If event_initial is set to true (the default is false) event iteration also occurs at the initial point, and at any point when resuming integration after a halt action has been performed.

The way in which sequences of events run is very relevant, and quite simply, they run in the order that they appear in the input events specification. If any event has triggered during an iteration, or the value of variables has changed (which criteria is in use is controlled by the event_iterate option), then the events are iterated through again, from the beginning, until event_maxiter is reached. The event_maxiter option is provided to prevent infinite event loops (which are surprisingly easy to create), and can be set to any positive integer value, but cannot be disabled.

There are two main ways in which event iteration can terminate. One of these is when a pass is made over all events, and none of the events trigger, and this is the default termination type (specified by event_iterate=trigger). The other termination type (event_iterate=prechange) compares the values of all variables (discrete and continuous) to their pre values (see below) after each complete pass over the events, and if no differences occur, the iteration terminates.

There are also several control actions that only make sense from within an event iteration:

The tobegin action exits the current event, and goes back to the start of the iteration (updating pre if relevant) while the toend action exits the current event, and jumps to the end of the iteration. Note that the only difference between these two behaviors is that when event_iterate=prechange, the toend action performs a termination check, while the tobegin action does not.

The nomark action is only useful for event_iterate=trigger, and exits the current event telling the control code to act as though the event did not execute, so as to not affect the trigger termination condition (e.g. if event_iterate=trigger, and the only action to run in an iteration ended with nomark, then the iteration would terminate).

The delayhalt action exits the current action and schedules a halt (a stop condition) for when the iteration terminates. This is in contrast to the halt action which exits the current action and halts immediately, exiting both the iteration and the continuous integration.

The breakiteration action exits the current event and immediately halts the event iteration, resuming continuous integration (unless a delayhalt was scheduled).

Another common mechanism required in event processing for simulation is the ability to examine a prior value of a variable. This can be used to check if a variable has changed, or a boolean variable has changed from true to false or false to true.

The mechanism available for this in dsolve[numeric] event processing is the pre operator. Note that pre has no meaning outside of event processing, as the value of pre is keyed by other events. The condition to detect if the value of, for example, u(t) has changed during event processing can be written as u(t) <> pre(u(t)).

There are some subtleties in dealing with pre, and these have to do with the point in the computation that the pre value refers to. There are four options for this, controlled by the event_pre optional argument:

iteration (default): pre(u(t)) refers to the value of u(t) at the start of the current iteration.

initial: pre(u(t)) refers to the value of u(t) at the start of the events at the current point (this value will never change during event iteration).

assignment: pre(u(t)) refers to the value of u(t) prior to the last assignment to u(t) since the start of the event iteration at the current point.

change: pre(u(t)) refers to the value of u(t) prior to the last nontrivial assignment to u(t) since the start of the events at the current point.

There are several simple aliases that are related to pre variables for use in the conditional part of the trigger of an event, and these are as follows:

change(u(t)): returns true if u(t) <> pre(u(t)), false otherwise.

risingedge(u(t)): returns true if pre(u(t)) is false and u(t) is true, or if pre(u(t)) < u(t) for non-boolean u(t).

fallingedge(u(t)): returns true if pre(u(t)) is true and u(t) is false, or if pre(u(t)) > u(t) for non-boolean u(t).

edge(u(t)): same as risingedge(u(t)).

In addition to the typical description of an event [trigger, action], additional options can be specified to control the behavior of events. Two options are currently available:

In addition to these options, a root-finding trigger (given as an expression) can be restricted to roots having a positive rate of change with respect to the independent variable (increasing) or a negative rate of change with respect to the independent variable (decreasing) by wrapping the root-finding trigger expression in a function with the name increasing or decreasing respectively. For example, the trigger increasing(y(t)) will only fire when y(t) passes through zero at a point t0 with D(y)(t0)>0.

There are many features available for dealing with events for procedure-form numerical ODE/DAE solutions that are an extension of the features previously available for stop conditions. This includes the abilities to query which event triggered a stop, continue numerical integration by clearing a continuous stop, and enable/disable events. The first example below demonstrates some of these features.

The features previously available through stop conditions include:

eventstop: query the last event that fired (0 if no event has fired)

eventclear: clear the last event allowing integration to continue

eventcount: the total number of events that have fired since integration started at the initial point

eventdisable=set: set of events to be disabled for further solution

eventenable=set: set of events to be enabled for further solution

eventstatus: query enabled/disabled events

In addition to the existing stop condition features a new event control and query feature, eventfired, has been added.

A call to a procedure-form dsolve numeric procedure with the argument eventfired=[n1,n2,n3], where n1,n2,n3 are numerical references to the location of a root-finding or time-interval event in the initial events specification list, will return a list [v1,v2,v3] of the times at which the event most recently fired. If the event has never fired, a time on the opposite side of the initial point (as the current integration direction) will be returned instead.

The calling sequence eventfired=[n1=v1,n2=v2,n3=v3] can be used to specify when each event has last fired (and as a result the events will not fire until one of the triggers fire pas the relevant specified value). When this form is used, the prior values of eventfired are returned in a list.

Any attempt to set or query eventfired for discrete or range conditions, or for an index outside of the bounds of the initial events list, will result in an error. Any attempt to set or query eventfired for a root-finding event with retrigger set to true will result in a warning as eventfired is not used for retriggerable events.

Until an integration direction is chosen (i.e. by requesting a solution value), eventfired cannot be used (as a point prior to the initial condition cannot be determined). It is important to note that the integration direction is an absolute direction relative to the initial value for the problem, not relative to the last computed point (dsolve procedures only ever integrate away from the initial condition for a problem). In order to accomplish this without requiring an initial computation, one can call the procedure with direction=1 to set a positive integration direction, or direction=-1 to set a negative integration direction. Note that setting the integration direction for a procedure that has already been used to compute a solution clears the prior computation, forcing recomputation from the initial point for any subsequent solution requests. If a direction cannot be determined (i.e. this is not done, and no computation has been performed) an error will result when calling with eventfired.

$\mathrm{dsn}≔\mathrm{dsolve}\left(\left\{\mathrm{diff}\left(y\left(t\right)\&comma;t\&comma;t\right)+y\left(t\right)=0\&comma;y\left(0\right)=0\&comma;\mathrm{D}\left(y\right)\left(0\right)=1\right\}\&comma;\mathrm{numeric}\&comma;\mathrm{events}=\left[\left[y\left(t\right)\&comma;\mathrm{halt}\right]\right]\right)$

$\mathrm{dsn}≔\mathbf{proc}\left(\mathrm{x\_rkf45}\right)...\mathbf{end\; proc}$

$\mathrm{dsn}\left(10\right)$

$\left[t=3.14159269672142\&comma;y\left(t\right)=1.73363927380432\times {10}^{\mathrm{-16}}\&comma;\frac{\&DifferentialD;}{\&DifferentialD;t}y\left(t\right)=\mathrm{-1.00000032537456}\right]$

$\mathrm{dsn}\left(\mathrm{eventclear}\right)$

$\mathrm{dsn}\left(10\right)$

$\left[t=6.28318539282259\&comma;y\left(t\right)=\mathrm{-3.99582710669533}\times {10}^{\mathrm{-16}}\&comma;\frac{\&DifferentialD;}{\&DifferentialD;t}y\left(t\right)=1.00000064490780\right]$

$\mathrm{dsn}\left(\mathrm{eventclear}\right)$

$\mathrm{dsn}\left(10\right)$

$\left[t=9.42477809037942\&comma;y\left(t\right)=\mathrm{-4.49049434789184}\times {10}^{\mathrm{-16}}\&comma;\frac{\&DifferentialD;}{\&DifferentialD;t}y\left(t\right)=\mathrm{-1.00000097453937}\right]$

$\mathrm{dsn}\left(\mathrm{eventclear}\right)$

$\mathrm{dsn}\left(10\right)$

$\left[t=10.\&comma;y\left(t\right)=\mathrm{-0.544021547556702}\&comma;\frac{\&DifferentialD;}{\&DifferentialD;t}y\left(t\right)=\mathrm{-0.839072448077328}\right]$

$\mathrm{dsn}\left(1.\times {10}^{\mathrm{-300}}\right)$

$\left[t=1.\times {10}^{\mathrm{-300}}\&comma;y\left(t\right)=1.00000000000000\times {10}^{\mathrm{-300}}\&comma;\frac{\&DifferentialD;}{\&DifferentialD;t}y\left(t\right)=1.\right]$

$\mathrm{dsn}\left(10\right)$

$\left[t=3.14159269672142\&comma;y\left(t\right)=1.73363927380432\times {10}^{\mathrm{-16}}\&comma;\frac{\&DifferentialD;}{\&DifferentialD;t}y\left(t\right)=\mathrm{-1.00000032537456}\right]$

$\mathrm{dsn}≔\mathrm{dsolve}\left(\left\{\mathrm{diff}\left(y\left(t\right)\&comma;t\&comma;t\right)+y\left(t\right)=0\&comma;y\left(0\right)=0\&comma;\mathrm{D}\left(y\right)\left(0\right)=1\right\}\&comma;\mathrm{numeric}\&comma;\mathrm{events}=\left[\left[\left[y\left(t\right)\&comma;0<\mathrm{diff}\left(y\left(t\right)\&comma;t\right)\right]\&comma;\mathrm{halt}\right]\&comma;\left[\left[y\left(t\right)\&comma;\mathrm{diff}\left(y\left(t\right)\&comma;t\right)<0\right]\&comma;\mathrm{halt}\right]\right]\right)$

$\mathrm{dsn}≔\mathbf{proc}\left(\mathrm{x\_rkf45}\right)...\mathbf{end\; proc}$

$\mathrm{dsn}\left(\mathrm{eventstatus}\right)$

$\mathrm{enabled}=\left\{1\&comma;2\right\},\mathrm{disabled}=\varnothing$

$\mathrm{dsn}\left(10\right)$

$\left[t=3.14159269672142\&comma;y\left(t\right)=1.73363927380432\times {10}^{\mathrm{-16}}\&comma;\frac{\&DifferentialD;}{\&DifferentialD;t}y\left(t\right)=\mathrm{-1.00000032537456}\right]$

$\mathrm{dsn}\left(\mathrm{eventstop}\right)$

$2$

$\mathrm{dsn}\left(\mathrm{eventclear}\right)$

$\mathrm{dsn}\left(10\right)$

$\left[t=6.28318539282259\&comma;y\left(t\right)=\mathrm{-3.99582710669533}\times {10}^{\mathrm{-16}}\&comma;\frac{\&DifferentialD;}{\&DifferentialD;t}y\left(t\right)=1.00000064490780\right]$

$\mathrm{dsn}\left(\mathrm{eventstop}\right)$

$1$

$\mathrm{dsn}\left(\mathrm{eventclear}\right)$

$\mathrm{dsn}\left(\mathrm{eventdisable}=\left\{2\right\}\right)$

$\mathrm{dsn}\left(10\right)$

$\left[t=10.\&comma;y\left(t\right)=\mathrm{-0.544021550832599}\&comma;\frac{\&DifferentialD;}{\&DifferentialD;t}y\left(t\right)=\mathrm{-0.839072452601256}\right]$

$\mathrm{dsn}\left(\mathrm{eventstop}\right)$

$0$

$\mathrm{dsn}\left(\mathrm{eventdisable}=\left\{1\right\}\right)$

$\mathrm{dsn}\left(\mathrm{eventstatus}\right)$

$\mathrm{enabled}=\varnothing ,\mathrm{disabled}=\left\{1\&comma;2\right\}$

$\mathrm{dsn}\left(100\right)$

$\left[t=100.\&comma;y\left(t\right)=\mathrm{-0.506372107554053}\&comma;\frac{\&DifferentialD;}{\&DifferentialD;t}y\left(t\right)=0.862327180873566\right]$

$\mathrm{dsn}\left(\mathrm{eventenable}=\left\{1\&comma;2\right\}\right)$

$\mathrm{dsn}\left(\mathrm{eventstatus}\right)$

$\mathrm{enabled}=\left\{1\&comma;2\right\},\mathrm{disabled}=\varnothing$

$\mathrm{dsn}\left(10\right)$

$\left[t=3.14159269672142\&comma;y\left(t\right)=1.73363927380432\times {10}^{\mathrm{-16}}\&comma;\frac{\&DifferentialD;}{\&DifferentialD;t}y\left(t\right)=\mathrm{-1.00000032537456}\right]$

$\mathrm{dsn}\left(\mathrm{eventstop}\right)$

$2$

$\mathrm{dsn}\left(-10\right)$

$\left[t=\mathrm{-3.14159269672142}\&comma;y\left(t\right)=\mathrm{-1.73363927380432}\times {10}^{\mathrm{-16}}\&comma;\frac{\&DifferentialD;}{\&DifferentialD;t}y\left(t\right)=\mathrm{-1.00000032537456}\right]$

$\mathrm{dsn}\left(\mathrm{eventstop}\right)$

$2$

$\mathrm{dsn}≔\mathrm{dsolve}\left(\left\{\mathrm{diff}\left(x\left(t\right)\&comma;t\&comma;t\right)=0\&comma;\mathrm{diff}\left(y\left(t\right)\&comma;t\&comma;t\right)=-9.8\&comma;x\left(0\right)=0\&comma;y\left(0\right)=1\&comma;\mathrm{D}\left(x\right)\left(0\right)=1\&comma;\mathrm{D}\left(y\right)\left(0\right)=0\right\}\&comma;\mathrm{numeric}\&comma;\mathrm{events}=\left[\left[y\left(t\right)\&comma;\mathrm{diff}\left(y\left(t\right)\&comma;t\right)=-\mathrm{diff}\left(y\left(t\right)\&comma;t\right)\right]\right]\&comma;\mathrm{range}=0..3\right)$

$\mathrm{dsn}≔\mathbf{proc}\left(\mathrm{x\_rkf45}\right)...\mathbf{end\; proc}$

$\mathrm{plots}\left[\mathrm{odeplot}\right]\left(\mathrm{dsn}\&comma;\left[x\left(t\right)\&comma;y\left(t\right)\right]\&comma;\mathrm{refine}=5\right)$

$\mathrm{dsys}≔\left\{\mathrm{diff}\left(\mathrm{x1}\left(t\right)\&comma;t\right)=1\&comma;\mathrm{diff}\left(\mathrm{x2}\left(t\right)\&comma;t\right)=-1\&comma;\mathrm{x1}\left(0\right)=-1\&comma;\mathrm{x2}\left(0\right)=1\right\}$

$\mathrm{dsys}≔\left\{\frac{\&DifferentialD;}{\&DifferentialD;t}\mathrm{x1}\left(t\right)=1\&comma;\frac{\&DifferentialD;}{\&DifferentialD;t}\mathrm{x2}\left(t\right)=\mathrm{-1}\&comma;\mathrm{x1}\left(0\right)=\mathrm{-1}\&comma;\mathrm{x2}\left(0\right)=1\right\}$

$\mathrm{dsn}≔\mathrm{dsolve}\left(\mathrm{dsys}\&comma;\mathrm{numeric}\&comma;\mathrm{events}=\left[\left[\mathrm{x1}\left(t\right)=1\&comma;\left[\mathrm{t1}=\mathrm{x1}\left(t\right)\&comma;\mathrm{x1}\left(t\right)=\mathrm{x2}\left(t\right)\&comma;\mathrm{x2}\left(t\right)=\mathrm{t1}\right]\right]\right]\&comma;\mathrm{range}=0..6\right)$

$\mathrm{dsn}≔\mathbf{proc}\left(\mathrm{x\_rkf45}\right)...\mathbf{end\; proc}$

$\mathrm{plots}\left[\mathrm{odeplot}\right]\left(\mathrm{dsn}\&comma;\left[t\&comma;\mathrm{x1}\left(t\right)\right]\&comma;\mathrm{refine}=1\right)$

$\mathrm{plots}\left[\mathrm{odeplot}\right]\left(\mathrm{dsn}\&comma;\left[t\&comma;\mathrm{x2}\left(t\right)\right]\&comma;\mathrm{refine}=1\right)$