Introduction

Differential equations alone are very effective for modeling continuous behavior of systems. However, many real systems involve components that change at discrete times, possibly triggered by states of the continuous solutions. For example, in a heating system, a thermostat will switch on once the temperature reaches a certain level.

An event trigger for a differential or differential algebraic system

is a point along the solution at which a Boolean-valued event function becomes True:

When is True for , then the value of can be found numerically during integration by evaluating with an approximation of the solution and using bisection.

In general, much more robust and faster ways of numerically finding the value are available if is of the form , where is a real-valued function. Symbolic processing is automatically done on the event expression to try to get a real-valued function or simpler form if possible. For example, the Boolean event is specially processed so that an event will trigger at time , and the Boolean event is automatically processed to the real function .

For some event actions, such as discontinuity crossings, it is necessary to find a bracketing interval such that is False and is True, or for a real-valued function , the signs of and are opposite.

When an event trigger is found, typically some action is performed. In some cases, the action can affect the integration of the system, by stopping integration or even changing the values of the solution variables, so the use of events in modeling gives a great deal of flexibility for discrete actions.

In NDSolve, events augment the differential equations with WhenEvent expressions.

This section shows basic examples using WhenEvent in NDSolve. Subsequent sections in the tutorial go into more depth on topics including event detection and location, actions that modify state variables, discontinuity handling, and application examples.

A simple example is locating an event, such as when a pendulum started at a nonequilibrium position will swing through its lowest point, and stopping the integration at that point.

From the solution you can see that the pendulum reaches its lowest point at about . Using the InterpolatingFunctionAnatomy package, it is possible to extract the value from the InterpolatingFunction object.

Event Detection and Location

Event location works in effectively two stages. In the event detection stage, for each time step taken by the numerical integration method each event is tested to see if it may have changed during the step. Once a possible event is detected, the location stage typically uses a root-finding procedure to find the actual event time .

Event Detection Methods

The simplest and fastest method for detecting whether an event occurred in a time step is to check for a change in the event function from False to True, or if it can be expressed as the root of a real function, a sign change.

The problem with just using the sign is that if a method takes large steps (such as extrapolation methods tend to) or the event function varies rapidly compared to the solution of the ODE, it is quite possible to miss events when double (or more) roots occur within a time step.

For events that can be expressed as a root of a real function, the method of using signs can be enhanced by including the time derivatives of the event function at the beginning and end of the step. For an ODE with event function , the time derivative is computed with

The values can be interpolated by a cubic Hermite polynomial that can be tested to see if it has any roots in the step interval. If necessary, the step is subdivided so that a root can be bracketed for the root-finding procedure.

Using the derivative makes the detection more robust with relatively little extra cost if does not depend on , so this is the default method for those cases. However, it is still entirely possible that events can be omitted using this method.

For some problems, it is vital to ensure that events are not missed, and NDSolve supports a method using interpolation that is very robust, but it can be significantly more expensive. This method is based on using the continuous or dense output from the method to create an interpolant of the event function with the same order as the method. The interpolant is tested for roots in the time step using interval arithmetic techniques.

The method for event detection is controlled by the "EventDetection" option of WhenEvent.

Settings for the "DetectionMethod" option.

When the event function varies much more rapidly than the solutions of the differential equation itself, the time steps may include many crossings of an event function. This can be mitigated by including the event function in the system as a dependent variable, using the equation for the time derivative of the event shown above. This can be done automatically using the option setting "IntegrateEvent"->True in WhenEvent.

When the event is integrated, the event interpolation is even more robust, since the event function interpolant is just a component of the dense output of the time integration method. In this case, the roots of the interpolant are the event locations to the full accuracy of the method, so there is no need for using a general root-finding method. For this reason, with "DetectionMethod"->"Interpolation", the event is integrated by default.

Event Location Methods

Once an event has been detected, the next task is to refine the location of the time within the time step where the event occurs. There are several different methods that can be used to refine the position. These include simply taking the solution at the beginning or the end of the integration interval, a linear interpolation of the event value, and using bracketed root-finding methods. The appropriate method to use depends on a trade-off between execution speed and location accuracy.

The method used to refine the location is controlled by the "LocationMethod" option of WhenEvent.

Settings for the "LocationMethod" option.

Location of a single event is usually fast enough so that the method used will not significantly influence the overall computation time. However, when an event is detected multiple times, the location refinement method can have a substantial effect.

"StepBegin" and "StepEnd" Methods

The crudest methods are appropriate for when the exact position of the event location does not really matter or does not reflect anything with precision in the underlying calculation.

Here is an example where "StepEnd" is appropriate because the location of the event is arbitrary.

Trigger an event to stop integration when a button is clicked. If the Stop button is clicked soon enough, the integration will be terminated before :

In this example, time steps are computed so quickly that there is no way that you can time the click of a button to be within a particular time step, much less at a particular point within a time step. Thus, based on the inherent accuracy of the event, there is no point in refining at all, so it is most appropriate to use the "StepBegin" or "StepEnd" location methods. In any example where the definition of the event is heuristic or somewhat imprecise, this is the best appropriate choice.

"LinearInterpolation" Method

When event results are needed for the purpose of points to plot in a graph, you only need to locate the event to the resolution of the graph. While just using the step end is usually too crude for this, a single linear interpolation based on the event function values suffices.

Denote the event function values at successive mesh points of the numerical integration:

Linear interpolation gives

A linear approximation of the event time is then:

Linear interpolation could also be used to approximate the solution at the event time. However, since derivative values and are available at the mesh points, a better approximation of the solution at the event can be computed cheaply using cubic Hermite interpolation as

for suitably defined interpolation weights:

You can specify refinement based on a single linear interpolation with the setting "LinearInterpolation".

At the resolution of the plot over the entire period, you cannot see that the endpoint may not be exactly where the derivative hits the axis. However, if you zoom in enough, you can see the error.

The linear interpolation method is sufficient for most viewing purposes, such as computing Poincaré sections. Note that for Boolean-valued event functions, linear interpolation is effectively only one bisection step, so the linear interpolation method may be inadequate for graphics.

Brent's Method

The default location method is the event location method "Brent", finding the location of the event using FindRoot as described in "Brent's Method". Brent's method starts with a bracketed root and combines steps based on interpolation and bisection, guaranteeing a convergence rate at least as good as bisection and in most cases much better. You can control the accuracy and precision to which FindRoot tries to get the root of the event function using method options for the "Brent" event location method. The default is to find the root to the same accuracy and precision as NDSolve is using for local error control.

When the event depends on the solutions of the differential equation, the argument needs to be approximated for within the time step. For methods that support continuous or dense output, the argument for the event function can be found quite efficiently simply by using the continuous output formula. However, for methods that do not support continuous output, the solution needs to be computed by taking a step of the underlying method, which can be relatively expensive. An alternate way of getting a solution approximation that is not accurate to the method order, but is consistent with using FindRoot on the InterpolatingFunction object returned from NDSolve, is to use cubic Hermite interpolation, obtaining approximate solution values in the middle of the step by interpolation based on the solution values and solution derivative values at the step ends.

For events that are set up to handle discontinuity crossings, it is very useful for the event handler to have a bracketing interval so that it can do proper evaluations on either side of the discontinuity and ensure the crossing is done correctly. This is done with the method "BrentBracketRoot".

Method options can be passed to the "Brent" and "BrentBracketRoot" methods that control the solution approximation and the use of FindRoot.

Options for event location methods "Brent" and "BrentBracketRoot".

If the AccuracyGoal or PrecisionGoal settings are Automatic, the value passed to FindRoot is based on the local error setting for NDSolve.

Location Method Comparison

This example integrates the pendulum equation for a number of different event location methods and compares the time when the event is found.

Event Actions

With WhenEvent[event,action], once an event trigger point is detected and located along a solution trajectory, then the event handler evaluates the action given in the second argument of WhenEvent.

First, the solution values are approximated, typically using a dense output formula for the time integration method. The derivatives are found by evaluating the right-hand side function for ODEs or using the dense output formula for DAEs. Then (in a way similar to Block) the time-independent variable is set to and the dependent variables and their derivatives are set to the approximate values in and , and the expression is evaluated. If the value returned from the evaluation is a Rule or a string directive, then the event handler will execute the directed action.

"StopIntegration"	stop integrating at ; return the solution
"RestartIntegration"	restart the integration at
"CrossDiscontinuity"	integrate to , extrapolate to , and restart at
"CrossSlidingDiscontinuity"	integrate to , extrapolate and check sliding mode conditions, and restart at
"RemoveEvent"	remove the event
y[t]->val	change the state variable y to val
d[t]->"DiscontinuitySignature"	change the discontinuity signature variable d

Event action directives.

With WhenEvent[event,{action₁,action₂,…}], action_i is evaluated and handled completely before action_i+1 is evaluated, so it is possible to give a sequence of actions that will be triggered for a single event.

There are many examples of event actions with and without directives in the reference page for WhenEvent. A few examples are given here that illustrate the sequential evaluation and the handling of directives.

Collect points where the and stop integrating the first time this occurs with <1:

Changing State Variables

In general, evaluation of an arbitrary action expression does not affect the solution state variables that are used by NDSolve. Using an action in WhenEvent that evaluates to a rule provides a flexible way of making state changes at a discrete event time.

A simple example of this is a model of a ball bouncing elastically on a flat surface.

The rule x'[t]->-r x[t] specifies that the state variable x^′[t] is replaced with -r x^′[t] whenever x[t]0. After the change is made, the integration is restarted, since effectively a new solution trajectory is being computed.

A change for several state variables can be specified either using a single rule with a list of variables and values or multiple rules. There are subtle semantic differences between the two approaches, demonstrated by the examples below.

In the rule {x[t],y[t]}->{y[t],x[t]}, the right-hand side is evaluated and then the corresponding state variables are set. When the action is a list of rules, each rule is evaluated in turn.

With the two rules, first is set to and then is set to the new value of , so the second rule actually has the effect of setting to .

For an ODE system , it is not possible to set the derivatives x'[t], since those are determined explicitly from the function. For higher-order systems that are automatically reduced to first order, this means that the highest-order derivative that appears in the equations cannot be set.

When solving as a system of differential-algebraic equations , it is possible to set the derivatives. A part of the DAE solving process is to find consistent values and so that the residual is 0. NDSolve is directed to treat equations as a system of DAEs if you use the Method option setting Method->{"EquationSimplification"->"Residual"}.

When either the variables or their derivatives are changed when solving as a DAE of the form , unless the new values given satisfy , a reinitialization needs to be computed where an attempt is made to find values of the variables that were not changed, so that . Details of the procedure and examples are given in "DAE Reinitialization with Events".

When the action is not explicitly a rule, the action is evaluated after replacing all explicit instances of Rule in the action expression with a special private head that has the attribute HoldFirst, so that the states on the left-hand side do not evaluate. If the return value has this special head, then a state change is executed.

At regular intervals, make a state change with a probability and show several random realizations of the simulation:

Discrete variables specified with the DiscreteVariables option can only be changed through rules. An advantage of using discrete variables is that they do not increase the system size that the differential equation solver methods work with—this is particularly important for stiff methods where the Jacobian is used.

Events and Discontinuous Differential Equations

The methods built into NDSolve for solving systems of differential equations are all effectively based on an implicit assumption that the differential equation satisfies some basic continuity conditions. For ODEs, a sufficient condition is Lipschitz continuity.

When a differential equation has discontinuities, these assumptions are violated and solvers may get poor quality solutions or, in some cases, may not be able to find solutions at all. Typically what happens is that in an attempt to satisfy local error tolerances, the solvers decrease step size until the step is so small that error tolerances appear to be satisfied even with the discontinuous function values, and the solver continues beyond the discontinuity with the approximation it found.

As an example, consider the ODE

where is a parameter.

Define a function for the right-hand side of the differential equation that switches across a discontinuity:

Solve the equation to for using the hidden definition counting the number of steps used:

The solution shown here was computed without any a priori knowledge of the discontinuity, and the reduction to the time step just prior to the discontinuity is seen. Once the discontinuity has been crossed, the time steps increase in size again. Essentially the same thing happens with the default (multistep) method, but a higher-order Runge–Kutta method was chosen for this example to emphasize the effect.

If the time at which a discontinuity is encountered can be identified, then an alternative way of handling the discontinuity is to step up to the discontinuity using values from the approaching side and then restart the integration using values on the other side of the discontinuity. Both the extremely small time steps and possible errors can be avoided this way.

A natural way of identifying the time at which a discontinuity is encountered is to set up an event. In the Wolfram Language, this can be done with WhenEvent. One issue with using events to detect discontinuities is that for events that depend on the solution, the solution may need to be computed past the event to correctly locate the event, but that involves the approximation of the solution across the discontinuity the event was set up to avoid in the first place. NDSolve avoids this problem by using discontinuity signature discrete variables that are only switched when the discontinuity is encountered, so that the solver in effect samples a smooth function. This is done automatically by NDSolve when it symbolically detects discontinuous functions in the differential equations, but it is possible to set up WhenEvent expressions explicitly that do the same thing; details are described following.

When NDSolve can see the symbolic form of the function, it is able to automatically identify the discontinuity and handle it using an event. This approach requires far fewer time steps.

The error makes a jump of several orders of magnitude when the discontinuity is not handled using an event.

When the parameter is increased sufficiently, numerical methods without events are not generally able to continue past the discontinuity because when the discontinuity is first encountered, the solution on either side points back toward the discontinuity.

Attempt to solve the equation to for using the hidden definition counting the number of steps used:

With , numerical methods without events are not generally able to continue past the discontinuity because when the discontinuity is first encountered, the solution on either side points back toward the discontinuity.

With the symbolic detection of the discontinuity, even though the Wolfram Language conditions for unique existence of the solution are violated, NDSolve is able to compute a continuation known as Filippov [F88] sliding mode in this situation. When the vector field points toward the discontinuity on both sides, the continuation effectively slides along the discontinuity surface until the field on one side or the other points away. This can be seen in the plot above where the solution is directly on the circle for some time, but eventually deviates.

The following sections of this tutorial explain the Filippov continuation and how to use discontinuity signatures.

Sliding Mode

When a discontinuity surface is encountered and the vector field on either side of the surface points back toward the solution, some sort of continuation method is necessary. Without a continuation, a numerical method will typically fail as it tries to take steps that cross and then recross the discontinuity; in some cases the numerical solution will display artificial chattering.

An appropriate continuation defined by Filippov [F88] effectively takes a convex combination of the components of the vector fields pointing along the discontinuity surface (or orthogonal to the gradient vector of a function for which the surface is a level curve).

Entering sliding mode.

In the preceding figure, the discontinuity surface is defined by and the vector fields for and are and , respectively. Then on the surface use , where is chosen so that is normal to the gradient , in particular . The solution can enter sliding mode either from above or below, as shown in the left and right parts of the figure, respectively.

The preceding definition is given for an autonomous system with no explicit dependence on the independent variable . The same ideas apply for a nonautonomous system by adding a variable and considering the autonomous system .

A simple example is provided by the vector field that has a discontinuity surface . For this vector field , , and ; therefore , so and the sliding motion is governed by .

The default method gets stuck because it eventually gets to tiny steps that chatter along the discontinuity.

The solution shown above is not correct; the larger oscillation or chattering of the solution near the sliding surface is an artifact of the numerical method.

When the discontinuous function is explicit in the equations, NDSolve can handle the discontinuity automatically by using the Filippov continuation.

The Filippov continuation is a natural choice since solving a sequence of systems with smooth functions approaching the discontinuous one approaches the sliding mode solution.

Define a smooth function that approaches as the parameter goes to zero:

Manipulate the solution over a range of :

In general, a trajectory will be trapped in sliding mode until one of the vector fields starts to point away. When this happens, the solution can be continued above or below the discontinuity.

Exiting sliding mode.

In the figure above, the discontinuity surface is defined by and the vector fields for and are and , respectively. Then on the surface use , where is chosen so that is normal to the gradient , in particular . Beyond the point shown in both sides of the figure, and both have the same sign. On the left because and , and on the right because and , and in either case the trajectory exits sliding mode. These points can be found by numerically finding either the root or as appropriate.

Consider the vector field . Once , the vector field shown following starts to point down again, and the trajectory can escape sliding mode.

It may occur that a solution enters and exits sliding mode repeatedly.

Consider the vector field . For , this is just an oscillator with negative damping. For , all trajectories approach the discontinuity surface . On , if is not too large, there is a set of points where there is sliding motion. The sliding motion in effect stabilizes the unbounded growth that would occur with the negative damping, and so there is a periodic orbit.

All of the examples shown so far are in two dimensions, since that is easier to visualize. However, the continuation works in any number of dimensions, and solutions can arise where the solution slides on more than one discontinuity surface at once.

Consider the two discontinuity surfaces in 3D defined by the surface of a sphere and the plane and the vector field given by .

The solution first slides along the surface of the sphere, then along both surfaces, and eventually just along the plane.

Discontinuity Signature

For functions with complicated definitions such that symbolic processing of discontinuities is not possible, it is possible to set up manually the same sort of discontinuity handling that NDSolve does automatically for discontinuous functions.

Discontinuity signature variables are set up by including a WhenEvent expression where the event is on the discontinuity surface, and the action is rule setting the variable to "DiscontinuitySignature" and using the DiscreteVariables option to specify the range of the variable to be either {-1,0,1} or {-1,1}.

Setting up a discontinuity signature variable s[t].

Consider the first example in the previous section on sliding mode with the vector field .

Plot the solution and the values of the discontinuity signature variable :

The variable switches from to 0 when the discontinuity is encountered and the solution goes into sliding mode. Sliding mode is indicated by a discontinuity signature of 0, since ideally the solution is exactly on the discontinuity surface and so is 0. In practice, the actual solution may not satisfy this due to numerical error. When in sliding mode, NDSolve tries to keep as close to the surface as possible.

If you know that a discontinuity will not lead to sliding mode, the needed computations can be done less expensively if you exclude 0 from the range of the discontinuity signature variable. However, if you do this and there could be sliding mode on the discontinuity, the solution will fail or be incorrect.

An autonomous system defined by

can be represented by , where is set by the Wolfram Language expression WhenEvent[e[x[t]]>0,s[t]->"DiscontinuitySignature"] in NDSolve.

The time derivative of the event function, is important in determining what happens when at some time . There are six cases where will be changed when with WhenEvent[e[x[t]]>0,s[t]->"DiscontinuitySignature"].

Crossing from below and above.

• Case 1: (Cross from below) and for with at . The solution trajectory simply crosses to and is set to .

• Case 2: (Cross from above) and for with at .The solution trajectory simply crosses to and is set to .



Entering sliding mode from above and below.

• Case 3: (Enter sliding mode from above) and for with at . Since is decreasing to zero it also holds that at . The derivative is negative above and positive below the discontinuity, so the trajectory is trapped on the discontinuity. This is referred to as sliding mode and the value of is set to 0. The effective vector field is computed using the Filippov continuation , where . Note that since and .

• Case 4: (Enter sliding mode from below) and for with at . This is essentially the same as case 3 and the solution enters sliding mode with set to 0.



Exiting sliding mode to above and below.

• Case 5: (Exit sliding mode to above) , , , and for and at with . The solution is no longer in sliding mode and is set to 1. Note that the event that triggers this change is effectively a zero crossing of .

• Case 6: (Exit sliding mode to below) , , , and for and at with . The solution is no longer in sliding mode and is set to . Note that the event that triggers this change is effectively a zero crossing of .

To see all of these cases in an actual example, consider the system of ODEs



where and are constants that can be defined to see the different cases.

Define a function that gives the differential equations given definitions fa for and fb for :
Define a function solver that, given function definitions and an initial condition for , solves the system using NDSolve with as a discontinuity signature variable:
Define a function that plots the solution superimposed on the vector field and a plot of the discontinuity signature variable along with if there is sliding mode:
Case 1: Cross from below with and :
Case 2: Cross from above with and :
Case 3: Enter sliding mode from above with and :
Case 4: Enter sliding mode from below with and :
Case 5: Exit sliding mode to above and :
Case 6: Exit sliding mode to below with and :

Examples

Falling Body

This system models a body falling under the force of gravity encountering air resistance (see [M04]).

Period of the Van der Pol Oscillator

The Van der Pol oscillator is an example of an extremely stiff system of ODEs. The event locator method can call any method for actually doing the integration of the ODE system. The default method, Automatic, automatically switches to a method appropriate for stiff systems when necessary, so that stiffness does not present a problem.

This integrates the Van der Pol system for a particular value of the parameter up to the point where the variable reaches its initial value and direction:

By selecting the endpoint of the NDSolve solution, it is possible to write a function that returns the period as a function of .

This defines a function that returns the period as a function of :

This uses the function to compute the period at :

This shows a log-log plot of the period as a function of :

Of course, it is easy to generalize this method to any system with periodic solutions.

Poincaré Sections

Using Poincaré sections is a useful technique for visualizing the solutions of high-dimensional differential systems.

For an interactive graphical interface see the package EquationTrekker.

The Hénon–Heiles System

Define the Hénon–Heiles system that models stellar motion in a galaxy.

The Poincaré section of interest in this case is the collection of points in the plane when the orbit passes through .

Since the actual result of the numerical integration is not required, it is possible to avoid storing all the data in InterpolatingFunction by specifying the output variables list (in the second argument to NDSolve) as empty, or {}. This means that NDSolve will produce no InterpolatingFunction as output, avoiding storing a lot of unnecessary data. NDSolve does give a message NDSolve::noout warning there will be no output functions, but it can safely be turned off in this case, since the data of interest is collected from the event actions.

The linear interpolation event location method is used because the purpose of the computation here is to view the results in a graph with relatively low resolution. If you were doing an example where you needed to zoom in on the graph to great detail or to find a feature, such as a fixed point of the Poincaré map, it would be more appropriate to use the default location method.

Specify a WhenEvent that defines the Poincaré section. The event action is to use Sow on the values of and :

This integrates the Hénon–Heiles system using a fourth-order explicit Runge–Kutta method with fixed step size of 0.25. The event action is to use Sow on the values of and :

Since the Hénon–Heiles system is Hamiltonian, a symplectic method gives much better qualitative results for this example.

This integrates the Hénon–Heiles system using a fourth-order symplectic partitioned Runge–Kutta method with fixed step size of 0.25. The event action is to use Sow on the values of and :

The ABC Flow

Bouncing Ball

This example is a generalization of an example in [SGT03]. It models a ball bouncing down a ramp with a given profile. The example is good for demonstrating how you can use multiple invocations of NDSolve with event location to model some behavior.

This defines the function for the bounce when the ball hits the ramp. The formula is based on reflection about the normal to the ramp, assuming only the fraction of energy is left after a bounce:

Avoiding Wraparound in PDEs

Many evolution equations model behavior on a spatial domain that is infinite or sufficiently large to make it impractical to discretize the entire domain without using specialized discretization methods. In practice, it is often the case that it is possible to use a smaller computational domain for as long as the solution of interest remains localized.

In situations where the boundaries of the computational domain are imposed by practical considerations rather than the actual model being studied, it is possible to pick boundary conditions appropriately. Using a pseudospectral method as described in "Pseudospectral Derivatives" with periodic boundary conditions can make it possible to increase the extent of the computational domain because of the superb resolution of the periodic pseudospectral approximation. The drawback of periodic boundary conditions is that signals that propagate past the boundary persist on the other side of the domain, affecting the solution through wraparound. It is possible to use an absorbing layer near the boundary to minimize these effects, but it is not always possible to completely eliminate them.

The sine-Gordon equation turns up in differential geometry and relativistic field theory. This example integrates the equation, starting with a localized initial condition that spreads out. The periodic pseudospectral method is used for the integration. Since no absorbing layer has been instituted near the boundaries, it is most appropriate to stop the integration once wraparound becomes significant. This condition is easily detected with event location using WhenEvent.

Three-Body Problem

This example illustrates the solution of the restricted three-body problem, a standard nonstiff test system of four equations. The example traces the path of a spaceship traveling around the Moon and returning to the Earth (see p. 246 of [SG75]). The ability to specify multiple events and the direction of the zero crossing is important.

The initial conditions have been chosen to make the orbit periodic. The value of corresponds to a spaceship traveling around the Moon and the Earth:

The first two solution components are coordinates of the body of infinitesimal mass, so plotting one against the other gives the orbit of the body.

Infection Model with Delay

Define a function that makes a WhenEvent expression to capture the local maximum of the ^th component. A separate tag for Sow is used to distinguish the components:

Discontinuity Signature and C Code

For symbolically specified differential equations, NDSolve will automatically scan the equations for discontinuous functions and set up events for handling the discontinuities appropriately. If a differential system is defined by evaluating a function that it is not possible to symbolically scan, then it may be necessary to set up discontinuity signature variables and events to handle discontinuities that may be embedded in the definition.

A function defined in C code that is compiled to link to the Wolfram Language with LibraryFunction is a good example of a function that the Wolfram Language cannot symbolically scan.

Solve the system for with the switching on the circle :

Show a plot of the solution with blue for , green for , and red for sliding mode: