NDSolve
NDSolve[eqns,u,{x,x_{min},x_{max}}]
finds a numerical solution to the ordinary differential equations eqns for the function u with the independent variable x in the range x_{min} to x_{max}.
NDSolve[eqns,u,{x,x_{min},x_{max}},{y,y_{min},y_{max}}]
solves the partial differential equations eqns over a rectangular region.
NDSolve[eqns,u,{x,y}∈Ω]
solves the partial differential equations eqns over the region Ω.
NDSolve[eqns,u,{t,t_{min},t_{max}},{x,y}∈Ω]
solves the timedependent partial differential equations eqns over the region Ω.
NDSolve[eqns,{u_{1},u_{2},…},…]
solves for the functions u_{i}.
Details and Options
 NDSolve gives results in terms of InterpolatingFunction objects.
 NDSolve[eqns,u[x],{x,x_{min},x_{max}}] gives solutions for u[x] rather than for the function u itself.
 Differential equations must be stated in terms of derivatives such as u'[x], obtained with D, not total derivatives obtained with Dt.
 Partial differential equations may also be specified using the differential operators Grad (∇), Div (∇.), Laplacian (∇^{2}), and Curl (∇). Typically these operators are used as in Inactive[op] to keep the operator form from evaluating.
 NDSolve solves a wide range of ordinary differential equations as well as many partial differential equations.
 NDSolve can also solve many delay differential equations.
 In ordinary differential equations, the functions u_{i} must depend only on the single variable t. In partial differential equations, they may depend on more than one variable.
 WhenEvent[event,action] may be included in the equations eqns to specify an action that occurs when event becomes True.
 The differential equations must contain enough initial or boundary conditions to determine the solutions for the u_{i} completely.
 Initial and boundary conditions are typically stated in the form u[x_{0}]==c_{0}, u'[x_{0}]==dc_{0}, etc., but may consist of more complicated equations.
 The c_{0}, dc_{0}, etc. can be lists, specifying that u[x] is a function with vector or general list values.
 Periodic boundary conditions can be specified using u[x_{0}]==u[x_{1}].
 The point x_{0} that appears in the initial or boundary conditions need not lie in the range x_{min} to x_{max} over which the solution is sought.
 Boundary values may also be specified using DirichletCondition and NeumannValue.
 In delay differential equations, initial history functions are given in the form u[x/;x<x_{0}]==c_{0}, where c_{0} is in general a function of x.
 The differential equations in NDSolve can involve complex numbers.
 NDSolve can solve many differential‐algebraic equations, in which some of the eqns are purely algebraic, or some of the variables are implicitly algebraic.
 The u_{i} can be functions of the dependent variables and need not include all such variables.
 The following options can be given:

AccuracyGoal Automatic digits of absolute accuracy sought Compiled Automatic whether expressions should be compiled automatically DependentVariables Automatic the list of all dependent variables EvaluationMonitor None expression to evaluate whenever the function is evaluated InitialSeeding {} seeding equations for some algorithms InterpolationOrder Automatic the continuity degree of the final output MaxStepFraction 1/10 maximum fraction of range to cover in each step MaxSteps Automatic maximum number of steps to take MaxStepSize Automatic maximum size of each step Method Automatic method to use NormFunction Automatic the norm to use for error estimation PrecisionGoal Automatic digits of precision sought StartingStepSize Automatic initial step size used StepMonitor None expression to evaluate when a step is taken WorkingPrecision MachinePrecision precision to use in internal computations  NDSolve adapts its step size so that the estimated error in the solution is just within the tolerances specified by PrecisionGoal and AccuracyGoal.
 The option NormFunction>f specifies that the estimated errors for each of the u_{i} should be combined using f[{e_{1},e_{2},…}].
 AccuracyGoal effectively specifies the absolute local error allowed at each step in finding a solution, while PrecisionGoal specifies the relative local error.
 If solutions must be followed accurately when their values are close to zero, AccuracyGoal should be set larger, or to Infinity.
 The default setting of Automatic for AccuracyGoal and PrecisionGoal is equivalent to WorkingPrecision/2.
 The default setting of Automatic for MaxSteps estimates the maximum number of steps to be taken by NDSolve, depending on start and stop time and an estimate of the step size. Should this not be possible, a fixed number of steps is taken.
 The setting for MaxStepFraction specifies the maximum step to be taken by NDSolve as a fraction of the range of values for each independent variable.
 With DependentVariables>Automatic, NDSolve attempts to determine the dependent variables by analyzing the equations given.
 NDSolve typically solves differential equations by going through several different stages, depending on the type of equations. With Method>{s_{1}>m_{1},s_{2}>m_{2},…}, stage s_{i} is handled by method m_{i}. The actual stages used and their order are determined by NDSolve, based on the problem to solve.
 Possible solution stages are:

"TimeIntegration" time integration for systems of differential equations "BoundaryValues" ordinary differential equation boundary value solutions "DiscontinuityProcessing" symbolic processing for handling of discontinuous differential equations "EquationSimplification" simplification of equation form for numerical evaluation "IndexReduction" symbolic index reduction for differential algebraic equations "DAEInitialization" consistent initialization for differential algebraic equations "PDEDiscretization" discretization for partial differential equations  With Method>m_{1} or Method>{m_{1},s_{2}>m_{2},…}, the method m_{1} is assumed to be for time integration, so Method>m_{1} is equivalent to Method>{"TimeIntegration">m_{1}}.
 Possible explicit time integration settings for the Method option include:

"Adams" predictor‐corrector Adams method with orders 1 through 12 "BDF" implicit backward differentiation formulas with orders 1 through 5 "ExplicitRungeKutta" adaptive embedded pairs of 2(1) through 9(8) Runge–Kutta methods "IDA" implicit backward differentiation formulas for DAEs "ImplicitRungeKutta" families of arbitrary‐order implicit Runge–Kutta methods "SymplecticPartitionedRungeKutta" interleaved Runge–Kutta methods for separable Hamiltonian systems  With Method>{"controller",Method>"submethod"} or Method>{"controller",Method>{m_{1},m_{2},…}}, possible controller methods include:

"Composition" compose a list of submethods "DoubleStep" adapt step size by the double‐step method "EventLocator" respond to specified events "Extrapolation" adapt order and step size using polynomial extrapolation "FixedStep" use a constant step size "OrthogonalProjection" project solutions to fulfill orthogonal constraints "Projection" project solutions to fulfill general constraints "Splitting" split equations and use different submethods "StiffnessSwitching" switch from explicit to implicit methods if stiffness is detected  Methods used mainly as submethods include:

"ExplicitEuler" forward Euler method "ExplicitMidpoint" midpoint rule method "ExplicitModifiedMidpoint" midpoint rule method with Gragg smoothing "LinearlyImplicitEuler" linearly implicit Euler method "LinearlyImplicitMidpoint" linearly implicit midpoint rule method "LinearlyImplicitModifiedMidpoint" linearly implicit Bader‐smoothed midpoint rule method "LocallyExact" numerical approximation to locally exact symbolic solution  The setting InterpolationOrder>All specifies that NDSolve should generate solutions that use interpolation of the same order as the underlying method used. »
Examples
open allclose allBasic Examples (7)
Solve a firstorder ordinary differential equation:
Use the function and its derivative in a plot:
Secondorder nonlinear ordinary differential equation:
Plot the function and its first two derivatives:
System of ordinary differential equations:
Solve the heat equation in one dimension:
Solve the Poisson equation over a Disk:
Find a minimal surface over a Disk with a sinusoidal boundary condition.
Solve a coupled nonlinear sineGordon equation over a region.
Scope (25)
Ordinary Differential Equations (8)
Specify any order equation. Reduction to normal form is done automatically:
Directly differentiate the solution to make a phase plot:
Directly specify a system of equations:
Solve for a vectorvalued function:
Plot the four components of the solution:
Different equivalent ways of specifying a harmonic oscillator as a secondorder equation:
As a system of firstorder equations:
Using a vector variable with the dimension deduced from the initial condition:
Use matrixvalued variables to compute the fundamental matrix solution:
Compare to the exact solution:
Define a Van der Pol equation:
The solution's "stiff" behavior that the default solver automatically handles:
Other methods may not be able to solve this system:
Equations can have multiple distinct solutions:
The solution y[x] is continuous, as it integrates the piecewise function once:
The solution y[x] is differentiable, whereas y'[x] is continuous only:
Partial Differential Equations (5)
Nonlinear advectiondiffusion equation in one dimension:
Define a system of PDEs of mixed parabolichyperbolic type:
Nonlinear sineGordon equation in two spatial dimensions with periodic boundary conditions:
Plot the solution at the final time:
Plot the time evolution of a radial cross section of the solution:
Solve a wave equation over a region with a slit:
Solve a Poisson equation with periodic boundary conditions on curved boundaries:
Boundary Value Problems (5)
A nonlinear multipoint boundary value problem:
Solve a nonlinear diffusion equation with Dirichlet and Neumann boundary conditions starting from an initial seed of .
Solve a nonlinear equation with Dirichlet boundary conditions starting from an initial seed of .
Solve a complexvalued nonlinear reaction equation with Dirichlet boundary conditions:
Delay Differential Equations (2)
Hybrid and Discontinuous Equations (4)
Options (29)
AccuracyGoal and PrecisionGoal (1)
DependentVariables (1)
EvaluationMonitor (2)
InitialSeedings (2)
InterpolationOrder (1)
Use InterpolationOrder>All to get interpolation the same order as the method:
This is more timeconsuming than the default interpolation order used:
MaxStepFraction (1)
Features with small relative size in the integration interval can be missed:
Use MaxStepFraction to ensure features are not missed, independent of interval size:
MaxSteps (1)
MaxStepSize (2)
The default step control may miss a suddenly varying feature:
A smaller MaxStepSize setting ensures that NDSolve catches the feature:
Attempting to compute the number of positive integers less than misses several events:
Setting a small enough MaxStepSize ensures that none of the events are missed:
Method (12)
TimeIntegration (4)
PDEDiscretization (2)
BoundaryValues (1)
DiscontinuityProcessing (1)
NDSolve automatically does processing for discontinuous functions like Sign:
If the processing is turned off, NDSolve may fail at the discontinuity point:
With some time integration methods, the solution may be very inaccurate:
An equivalent way to find the solution is to use "DiscontinuitySignature":
The solutions are effectively identical:
The discontinuity signature is 0 when the solution is in sliding mode:
EquationSimplification (2)
The solution cannot be completed because the square root function is not sufficiently smooth:
One solution can be found by forming a residual and solving as a DAE system:
The other solution branch can be given by specifying a consistent value of :
With the suboption "SimplifySystem">True, NDSolve uses symbolic solutions for components with a sufficiently simple form:
IndexReduction (1)
An index 3 formulation of a constrained pendulum using index reduction:
The default method can only solve index 1 problems:
The problem resulting from symbolic index reduction can be solved:
Solve using reduction to index 0 and a projection method to maintain the constraints:
Plot implicit energy constraint for the two solutions at the time steps:
DAEInitialization (1)
Use forward collocation for initialization to avoid problems with the Abs term at 0:
NormFunction (1)
StartingStepSize (1)
StepMonitor (3)
WorkingPrecision (1)
Error in the solution to a harmonic oscillator over 100 periods:
When the working precision is increased, the local tolerances are correspondingly increased:
With a large working precision, sometimes the "Extrapolation" method is quite effective:
Error in the solution to a harmonic oscillator over 100 periods:
When the working precision is increased, the local tolerances are correspondingly increased:
With a large working precision, sometimes the "Extrapolation" method is quite effective:
Applications (31)
Ordinary Differential Equations (5)
Simulate Duffing's equation for a particle in a double potential well:
The solution depends strongly on initial conditions:
The Lorenz equations [more info]:
The Lotka–Volterra predatorprey equations [more info]:
Look at the appearance of the blue sky catastrophe orbit in the Gavrilov–Shilnikov model:
Reduced 3body problem [more info]:
A formulation suitable for a number of different initial conditions:
Partial Differential Equations (6)
Simple model for soil temperature at depth x with periodic heating at the surface:
Simple wave evolution with periodic boundary conditions:
Wolfram's nonlinear wave equation [more info]:
Wolfram's nonlinear wave equation in two space dimensions:
A soliton profile perturbed by a periodic potential in a nonlinear Schrödinger equation:
Use Stokes's equation to compute the fluid velocity field in a narrowing channel:
Delay Differential Equations (1)
Hybrid Differential Equations (5)
Simulate a bouncing ball that retains 95% of its velocity in each bounce:
Model a ball bouncing down steps:
Each time a linear oscillator solution crosses the negative axis, reflect it across the axis:
The solution of this reset oscillator exhibits chaotic behavior:
Plot the solution on the negative axis with a histogram of the reflection points:
Model a onedegreeoffreedom impact oscillator with sinusoidal forcing:
Model a damped oscillator that gets a kick at regular time intervals:
Mechanical Systems (3)
Model the motion of a pendulum in Cartesian coordinates. Derive the governing equations using Newton's second law of motion, and , with a force diagram:
Add damping to the pendulum so that it slows down over time:
The pendulum stabilizes at the vertical fixed point:
Change the rod to a stiff spring by modifying the constraint :
The solution contains highfrequency spring oscillations:
Model a double pendulum of unit mass and length that is released from the horizontal plane. Begin by deriving the equations of motion using Newton's second law:
Simulate the system by enforcing the equations and constraints as invariants:
Visualize the motion of the double pendulum:
Model a block on a moving conveyor belt anchored to a wall by a spring using different models for the friction force , including viscous, Coulomb, Stribeck, and static. Compare positions and velocities for the different models:
Newton's equation for the block:
Viscous friction is proportional to the relative velocity :
The block stabilizes just above the spring's natural length of 1:
Coulomb friction is proportional to the sign of relative velocity :
The block moves with the belt until the spring force is strong enough:
Stribeck friction is a refined Coulomb friction :
The variation at low velocities is slightly reduced:
Static friction holds the block in place until the spring force exceeds some value μ depending on roughness of surfaces. Use the discrete variable stuck set to 1 when the block is stuck and 0 otherwise:
Check whether the spring force is smaller than μ, and if the block is not moving relative to the belt:
The block repeatedly sticks to the belt, then slips away due to the spring force:
Electrical Systems (6)
Simulate the response of an RLC circuit to a step in the voltage at time :
Use component laws together with Kirchhoff's laws for connections:
Simulate the response of an RLC circuit to a step in the voltage at time :
Simulate the behavior of a parallel RLC circuit:
Show the response under a constant input current:
Show the currents in the R, L, C components and the resulting voltage:
Model a transistoramplifier circuit:
The input voltage varies sinusoidally:
The transistor dispatches the voltage in a nonlinear way, depending on :
Use Ohm's law and Kirchhoff's law to determine the governing equations for each node:
Simulate the singular system of equations:
The transistor amplifies voltage relative to :
Model a DCtoDC boost converter from input voltage level vi to desired output voltage level vd using a pulsewidth modulated feedback control q[t]:
Use Kirchhoff's laws to get a model for the circuit above:
The control signal q[t] will switch the transistor on for a fraction vi/vd of each period τ:
Boost from a lower voltage vi=24 to a higher voltage vd=36:
Model a DCtoDC buckboost converter from input voltage level vi to desired output voltage level vo using a pulsewidth modulated feedback control q[t]:
Use Kirchhoff's laws to get a model for the circuit above:
The control signal q[t] will switch the transistor on for a fraction vd/(vi+vd) of each period:
Hydraulic Systems (3)
Model the change in height of water in two cylindrical tanks as water flows from one tank to another through a pipe:
Use pressure relations and mass conservation:
Model the flow across the pipe with the Hagen–Poiseuille relation:
Model the second tank as a leaky tank:
Due to the leak in the second tank, both tanks will eventually drain out:
Model the change in height of water in two hemispherical tanks as water flows from one tank to another through a pipe:
Model the change in height of water in three tanks such that one tank feeds water to the other tanks:
Model the flow across the pipe with the Hagen–Poiseuille relation:
The flow rate from the first pipe will equal the sum of the flow rates in other two pipes:
Model the third tank as a leaky tank:
The first two tanks reach equilibrium and then drain at the same rate:
Chemical Systems (2)
Model the kinetics of an autocatalytic reaction:
The rate equations are given as:
The concentration of the species a, b, c should always be a constant:
Solve and visualize the evolution of the three species:
Model a chemical process of two species, FLB and ZHU, that are continuously mixed with carbon dioxide:
The inflow of carbon dioxide per volume unit is denoted by:
The rate equations are given as:
Equilibrium equation is given as:
Solve the equations and determine the concentration change in FLB, ZHU, CO2, and ZLA:
Properties & Relations (7)
Symbolic versus numerical differential equation solving:
The defining equation for JacobiSN:
Numerically compute values of an integral at different points in an interval:
For functions of the independent variable, NDSolve effectively gives an indefinite integral:
Finding an event is related to finding a root of a function of the solution:
Event location finds the root accurately and efficiently:
This gives as a function of for a differential equation:
Solve the equivalent boundary value problem:
Use NDSolve as a solver for a SystemModel:
Plot variables from the simulation result:
Use SystemModel to model larger hierarchical models:
Possible Issues (11)
Numerical Error (3)
The error tends to grow as you go further from the initial conditions:
Find the difference between numerical and exact solutions:
Error for a nonlinear equation:
For highorder methods, the default interpolation may have large errors between steps:
Interpolation with the order corresponding to the method reduces the error between steps:
Differential Algebraic Equations (3)
NDSolve cannot automatically handle systems of index greater than 1:
Highindex systems can be solved by performing index reduction on the system:
Here is a system of differentialalgebraic equations:
NDSolve may change the specified initial conditions if it cannot find the solution with :
Change the initial starting guess for the iterations to avoid such issues:
NDSolve is limited to index 1, but the solution with has index 2:
To solve highindex systems, use index reduction to reduce the DAE to index 1:
The default method may not be able to converge to the default tolerances:
With lower AccuracyGoal and PrecisionGoal settings, a solution is found:
The "StateSpace" time integration method can solve this with default tolerances:
Partial Differential Equations (3)
The spatial discretization is based on the initial value, which varies less than the final value:
By increasing the minimal number of spatial grid points, you can accurately compute the final value:
The plot demonstrates the onset of a spatially more complex solution:
Define a heat equation with an initial value that is a step function:
Discontinuities in the initial value may result in too many spatial grid points:
Setting the number of spatial grid points smaller results in essentially as good a solution:
Define a Laplace equation with initial values:
The solver only works for equations well posed as initial value (Cauchy) problems: