The dsolve command with the numeric or type=numeric option on a real-valued two-point boundary value problem (BVP), finds a numerical solution for the ODE or ODE system BVP.

The type of problem (BVP or IVP) is automatically detected by dsolve, and an applicable algorithm is used. The optional equation method=bvp[submethod] indicates that a specific BVP solver is to be used. The available submethods are traprich, trapdefer, midrich, and middefer. The first two methods, traprich and trapdefer,are trapezoid methods that use Richardson extrapolation enhancement or deferred correction enhancement, respectively. The remaining two methods, midrich and middefer, are midpoint methods with the same enhancement schemes.

There are two major considerations when choosing a submethod for a BVP. The trapezoid submethods are generally more efficient for typical problems, but the midpoint submethods are capable of handling harmless end-point singularities that the trapezoid submethods cannot. For the enhancement schemes, Richardson extrapolation is generally faster, but deferred corrections uses less memory on difficult problems.

The available methods are fairly general, and should work on a variety of BVPs. They are not well suited to solve BVPs that are stiff or have solutions with singularities in their higher order derivatives.

The solution method is capable of handling both linear and nonlinear BVPs with fixed, periodic, and even nonlinear boundary conditions. For nonlinear boundary conditions, however, it may be necessary to provide an initial solution profile that approximately satisfies the boundary conditions (see 'approxsoln' below). The method is also capable of handling BVP systems with undetermined parameters, given a sufficient number of boundary conditions to determine their values.

Computation can be performed in both hardware precision and arbitrary precision, based on the setting of Digits. If Digits is smaller than the hardware precision for the machine, then computations are performed in hardware precision (see evalhf). If Digits is larger, then computations are performed in Maple floating point. In both cases, many of the more computationally intensive steps are performed in compiled external libraries.

The return value of dsolve and the manipulation of the input system is controlled by the following three optional equations, which are discussed in dsolve[numeric].

In addition to the output options available for all numerical ODE problems, the option output=mesh is also available. It specifies that dsolve return an array-form output where independent variable values are chosen as the discrete mesh used internally by the method. This also has the effect of changing the interpolant default to false because only the discrete solution is required. See interpolant below.

The 'known' option specifies user-defined known functions, and is discussed in dsolve[numeric].

The options listed below are specific to BVP and include options to control the solution process, the accuracy of the result, and the initial mesh or solution profile to be used.

'abserr'= numeric

Numeric value that gives an absolute error tolerance for the solution. In all but exceptional cases the true solution should be within the error tolerance value of the continuous approximate solution obtained by the method.  The default value of abserr is Float(1,-6).

Note: This is an absolute error tolerance, so if the scale of the problem is such that the solution values are much smaller/larger than unitary, this parameter requires adjustment.

As for IVPs (see dsolve[numeric,IVP]), the system is converted to a first order system internally before the numerical solution is computed. The absolute error tolerance is applied to all components of the converted system.

'range'= numeric..numeric

Gives the left and right boundary points of the solution interval in the form leftpt..rightpt. This enables the method to be used to compute the solution of an initial value problem over a fixed interval using the absolute global error bound specified by abserr.

'adaptive'= boolean

Boolean value that determines whether mesh adaptation is used to obtain the solution. By default this is true.

The method applied is based on arc-length, with additional rules to prevent adjacent steps from changing their size too rapidly, restrictions on the largest and smallest allowable step sizes as compared with the corresponding fixed step size, and restrictions preventing the mesh points near the boundaries from becoming too widely spaced.

Note: Some problems work better with fixed step-size meshes.

'maxmesh'= integer

Integer value that determines the maximum number of points dsolve uses to compute the numerical solution. The numeric BVP solver internally uses a discrete mesh of points to calculate the approximate solution, which is adjusted as greater accuracy is required. If the desired accuracy cannot be obtained with the current limitation imposed by maxmesh, an error is returned. The default value for maxmesh is $128$. Its value must be between $32$ and $134217728$.

'initmesh'= integer

Integer value that determines the number of points dsolve uses to compute the initial solution profile. In some cases, the default initial $8$ point mesh does not have sufficient resolution to obtain the initial solution profile, so increasing this value can give a solution when the default value does not. Its value must be between $8$ and $134217728$.

'approxsoln'= array, list or procedure

Argument that enables specification of an initial approximate discrete solution to be used as a starting point for the problem. There are many forms that this argument can take.

The simplest forms are the output of another dsolve/numeric computation when output is specified as an array (Matrix form output), or as a procedure (procedurelist, listprocedure, or operator output).

Another form is as a two-dimensional array where the names of the dependent variables and independent variables are given in the first column, followed by their values at different mesh points in the following columns. All dependent variables in the (converted to first order) system, must be present. This form can also be specified using nested lists.

The final form is as a list of equations that describe the initial solution profile as functions of the independent variable. For example, for a second order BVP in $u\left(x\right)$, this could be specified as $\left[u\left(x\right)=f\left(x\right)\right]$ or as $\left[u\left(x\right)=f\left(x\right)\,\frac{\ⅆ}{\ⅆx}u\left(x\right)=g\left(x\right)\right]$, where the $f\left(x\right),g\left(x\right)$ are fully determined functions of $x$. In the former case, the required derivative values of $u\left(x\right)$ are computed by evaluation of the derivative of $f\left(x\right)$.

In the cases with a discrete solution given in array form, the values of the independent variable must always be provided at either boundary, in increasing order over the solution region, and at a minimum of $8$ distinct points. In addition, the specification of a discrete approximate solution cannot be used in combination with initmesh.

For systems with free parameters (determined through interaction with the boundary conditions) it is suggested that approximate values for the parameters also be provided.

'continuation'= name

Argument that allows solution of a BVP via a continuous transformation from an easier problem to the desired problem. This method is used in obtaining an initial solution profile, and is most helpful for problems for which the Newton iteration for the initial solution approximation does not converge. The continuation option provides the name for the continuation parameter for the problem. The continuation problem must be constructed so that the parameter c, when varied from $0$ to $1$, defines a different BVP for each value of c, where $c=0$ represents the simpler problem, and $c=1$ represents the problem to be solved. The continuation parameter can be present in the differential equation, the values for the boundary conditions, or both. Two examples of the use of continuation for the solution of BVP problems can be found in dsolve[numeric_bvp,advanced].

For best results, the continuation problem should be constructed so that the BVP solutions vary continuously with the parameter, and so that the rate of change of the solution, with respect to the parameter, is very roughly constant over the range $c=0..1$.

'mincont'= numeric

Argument that provides a minimum value for the allowed change in the continuation parameter when computing the initial solution profile. The value must be positive, and less than or equal to $\frac{1}{10}$. The default value is $\frac{1}{100}$. This option is valid only when continuation has also been specified.

For some discussion of difficult BVP problems, and possible solutions, see dsolve[numeric_bvp,advanced].

'interpolant'= boolean

Boolean that is true by default, and specifies that the solution should be refined until an interpolation of the solution gives the requested accuracy. When this is false, the solution is improved until only the discrete solution has the requested accuracy. Also see the output=mesh argument.

Results can be plotted by using the function odeplot in the plots package.

$\mathrm{dsol1}≔\mathrm{dsolve}\left(\left\{\mathrm{diff}\left(y\left(x\right)\,x\,x\right)-2y\left(x\right)=0\,y\left(0\right)=1.2\,y\left(1\right)=0.9\right\}\,\mathrm{numeric}\right)$

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

$\mathrm{dsol1}\left(0\right)$

$\left[x=0.\,y\left(x\right)=1.20000000000000\,\frac{\ⅆ}{\ⅆx}y\left(x\right)=\mathrm{-1.25251895272792}\right]$

$\mathrm{dsol1}\left(0.2\right)$

$\left[x=0.2\,y\left(x\right)=0.994463627112648\,\frac{\ⅆ}{\ⅆx}y\left(x\right)=\mathrm{-0.816528958014159}\right]$

$\mathrm{dsol1}\left(0.5\right)$

$\left[x=0.5\,y\left(x\right)=0.832942089176703\,\frac{\ⅆ}{\ⅆx}y\left(x\right)=\mathrm{-0.276385195570455}\right]$

$\mathrm{dsol1}\left(1\right)$

$\left[x=1.\,y\left(x\right)=0.900000000000000\,\frac{\ⅆ}{\ⅆx}y\left(x\right)=0.555701111490340\right]$

$\mathrm{deq2}≔\mathrm{diff}\left(y\left(x\right)\,x\,x\right)+\frac{\left(2+{y\left(x\right)}^2\right)y\left(x\right)}{1+{y\left(x\right)}^2}=1$

$\mathrm{deq2}≔\frac{{\ⅆ}^2}{\ⅆx^2}y\left(x\right)+\frac{\left(2+{y\left(x\right)}^2\right)y\left(x\right)}{1+{y\left(x\right)}^2}=1$

$\mathrm{bc2}≔y\left(0\right)=0,y\left(2\right)=3\:$

$\mathrm{dsol2}≔\mathrm{dsolve}\left(\left\{\mathrm{bc2}\,\mathrm{deq2}\right\}\,\mathrm{numeric}\,\mathrm{output}=\mathrm{listprocedure}\right)$

$\mathrm{dsol2}≔\left[x=\mathbf{proc}\left(x\right)...\mathbf{end\; proc}\,y\left(x\right)=\mathbf{proc}\left(x\right)...\mathbf{end\; proc}\,\frac{\ⅆ}{\ⅆx}y\left(x\right)=\mathbf{proc}\left(x\right)...\mathbf{end\; proc}\right]$

$\mathrm{fy}≔\mathrm{subs}\left(\mathrm{dsol2}\,y\left(x\right)\right)$

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

$\left[\mathrm{seq}\left(\mathrm{fy}\left(\frac{2i}{4}\right)\,i=0..4\right)\right]$

$\left[0.\,1.20426489602605\,2.24834461794798\,2.89441120924878\,3.00000000000000\right]$

$\mathrm{dsys3}≔\mathrm{diff}\left(y\left(x\right)\,x\right)-y\left(x\right)=0,y\left(1\right)=1$

$\mathrm{dsys3}≔\frac{\ⅆ}{\ⅆx}y\left(x\right)-y\left(x\right)=0,y\left(1\right)=1$

$\mathrm{dsol3}≔\mathrm{dsolve}\left(\left\{\mathrm{dsys3}\right\}\,\mathrm{numeric}\,\mathrm{method}=\mathrm{bvp}\,\mathrm{range}=0..1\,\mathrm{abserr}=1.\times {10}^{\mathrm{-10}}\right)$

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

$\mathrm{dsol3}\left(0\right)$

$\left[x=0.\,y\left(x\right)=0.367879441157275\right]$

$\mathrm{evalf}\left[11\right]\left(\exp \left(-1\right)\right)$

$0.36787944117$

$\mathrm{dsys4}≔\mathrm{diff}\left(y\left(x\right)\,x\,x\right)-ay\left(x\right)=0,y\left(0\right)=1,y\left(1\right)=1,\mathrm{D}\left(y\right)\left(0\right)=2$

$\mathrm{dsys4}≔\frac{{\ⅆ}^2}{\ⅆx^2}y\left(x\right)-ay\left(x\right)=0,y\left(0\right)=1,y\left(1\right)=1,\mathrm{D}\left(y\right)\left(0\right)=2$

$\mathrm{dsol4}≔\mathrm{dsolve}\left(\left\{\mathrm{dsys4}\right\}\,\mathrm{numeric}\,\mathrm{output}=\mathrm{operator}\right)$

$\mathrm{dsol4}≔\left[x=\mathbf{proc}\left(x\right)...\mathbf{end\; proc}\,y=\mathbf{proc}\left(x\right)...\mathbf{end\; proc}\,\mathrm{D}\left(y\right)=\mathbf{proc}\left(x\right)...\mathbf{end\; proc}\,a=\mathbf{proc}\left(x\right)...\mathbf{end\; proc}\right]$

$\mathrm{dsol4}\left(0\right)$

$\left[x=0\,y\left(0\right)=1.00000000000000\,\mathrm{D}\left(y\right)\left(0\right)=2.00000000000000\,a\left(0\right)=\mathrm{-2.96069553617247}\right]$

$\mathrm{dsol4}\left(0.5\right)$

$\left[x=0.5\,y\left(0.5\right)=1.53330815132594\,\mathrm{D}\left(y\right)\left(0.5\right)=5.80015080018856\times {10}^{\mathrm{-11}}\,a\left(0.5\right)=\mathrm{-2.96069553617247}\right]$

Ascher, U.; Mattheij, R.; and Russell, R. "Numerical Solution of Boundary Value Problems for Ordinary Differential Equations." SIAM Classics in Applied Mathematics. Vol. 13. (1995).

Ascher, U., and Petzold, L. "Computer Methods for Ordinary Differential Equations and Differential-Algebraic Equations." SIAM, Philadelphia. 1998.