The pdsolve(PDEsys,conditions,numeric,vars,options) command returns a module that can be used to compute numerical solutions for time-based PDE systems over a fixed finite 1-space interval.

The pdsolve/numeric routine uses finite difference methods to obtain these numerical solutions.

The PDE system PDEsys must be specified in diff notation.

Derivatives in the initial/boundary conditions conditions must be specified in indexed D notation (for example, for u(x,y), ${\mathrm{D}}_1\left(u\right)\left(0\,y\right)$ describes the normal derivative on the line x=0).

The vars parameter specifies the dependent variable(s) of the problem. For the solvers that can compute solutions of systems of PDEs, it provides an order for the dependent variables in the output.

For information on specifying the independent variables of the problem, see the Options section of this help page.

The solver has two modes of operation.

- The first mode of operation uses the default method, which is a centered implicit scheme, and is capable of finding solutions for higher order PDE or PDE systems. The PDE systems must be sufficiently close to a standard form for the method to find the numerical solution.

- The second mode of operation is strictly educational, and allows specification of a particular method to solve a PDE. This mode is restricted to a single PDE.

Note: This page describes only the default method. For information on the specified method approach, see Finding a Numerical Solution of a PDE Using a Specified Method. Both modes accept Dirichlet, Neumann, Robin, or periodic boundary conditions.

The optional equations for the default method are as follows.

'indepvars'=list

Specifies the independent variables of the problem, and provides an order to the independent variables for returned procedures that compute the solution given the values of the independent variables.

'time'=name, 'range'=l..r

The time argument specifies the name of the time variable for the problem, while the range argument specifies the spatial domain of the problem. Note: The values of both options are determined automatically if the boundary conditions are specified for both end points of the domain. These options are typically needed only for first order non-periodic problems where only one of the end points appears in the condition.

'spacestep'=numeric

Specifies the spacing of the spatial points on the discrete mesh on which the solution is computed, and defaults to 1/20th of the spatial range of the problem, $r-l$, where $l$ is the left boundary, and $r$ is the right. Note: The spacing must be small enough that a sufficient number of points are in the spatial domain for the given method, boundary conditions, and spatial interpolation (see below). If the given value of spacestep does not fit into the input domain an integral number of times, the closest smaller step size that does is used. For problems using error estimates, or an adaptive approach, the total number of spatial points must be odd, so again, if this is not the case for the specified spacestep, then it is reduced to satisfy this requirement.

'timestep'=numeric

Provides the size of the time step used when integrating the PDE in the time direction, and defaults to the value of the spatial step spacestep before adjustments are made for fitting the interval.

'bcpts'=integer

Specifies the number of points used in discretization of the boundary conditions of the problem. It defaults to one greater than the differential order of the PDE in the space variable. This provides a first order accurate approximation to the boundary condition. Additional information is available below in Details.

'optimize'=boolean/symbol

Specifies the level of optimization to use in generating the method. Three types of optimization are available, true, partial, and false. The default partial optimization gives fairly good performance for a low time penalty in generating the method. A setting of false disables optimization altogether (not recommended). A setting of true causes pdsolve to call codegen[optimize] when constructing the numerical method. Because the additional time efficiency from the use of the fully optimized method is often offset by the time taken to perform the optimization, partial is the default setting.

'errorest'=boolean,'errortype'=keyword,'abstol'=numeric,'mintimestep'=numeric, 'maxtimestep'=numeric

These are all options related to error estimation, and error control, and are discussed in the pdsolve,numeric,errorcontrol page.

'compile'=boolean

This option specifies that the underlying procedures used to evaluate the PDE should be compiled. By default this is false, as not all PDE can be compiled, and this option can only be used for the default PDE solution method. Problems typically run significantly faster when compiled. Note that in the event the compile fails, the problem will be silently run with compile=false.

The default method uses a second order (in space and time) centered, implicit finite difference scheme to obtain the computed solution. The number of points in the stencil of the finite difference scheme is one greater than the order of each equation.

Note: A first order accurate accurate boundary condition (the default) can be used for second order accurate method, and still provide a solution of second order accuracy. Use of boundary conditions with an order of accuracy that is too low degrades the accuracy of the solution, while use of boundary conditions with an order of accuracy that is too high may negatively impact the stability of the problem. For first and second order accurate methods, it is generally recommended to use one or two more boundary points than the differential order of the PDE in the space variable. For third and fourth order accurate methods it is generally recommended to use two or three more boundary points than the differential order of the PDE in the space variable.

PDE that are greater than first order in time are solved by converting to an equivalent first order system. Note: The conversion to first order is the expected one, namely that for any dependent variables that are higher than first order, variables representing all lower order time derivatives, and corresponding equations, are appended to the system. This choice of additional variables may not always be desirable, particularly when error control is in use, so often better performance can be obtained by choosing the first order system yourself. For more information on this topic, see pdsolve/numeric/firstorder.

Requested values that do not lie on a mesh point are computed using a Lagrange interpolation that uses 4 space points and 3 time points. The time points are chosen (whenever possible) so that only one time value is larger than the requested time.

The module returned by pdsolve/numeric has a number of methods that can be used to compute or view the solution of the input PDE: plot, plot3d, animate, and value. In addition, the settings method is available to query/set certain parameters of the solution process.

If the method returned by pdsolve/numeric is assigned to a name, say pde_module, the calling sequences that use the methods have the form pde_module:-method(arguments), where method is one of plot, plot3d, animate, value, and settings.

depspec

All methods involving the computation of solution values allow specification of multiple functions of the dependent and independent variables along with options specific to each function. This allows multiple plots to be displayed with the same command using different options for each. The full description of this argument is here. In the methods below, it is only referenced by name.

The default value of depspec is the first dependent variable of the problem with no additional options.

In the single function case, depspec can be a function of the dependent and independent variables, or a list consisting of a function of the dependent and independent variables followed by options. For example, $u\left(x\,t\right)-\sin \left(x\right)$, or $\left[u\left(x\,t\right)\,\mathrm{color}=\mathrm{blue}\right]$.

If specifying multiple functions, depspec can be given as a list of functions or function/option combinations as described in the single function case. For example, $\left\{u\left(x\,t\right)\,w\left(x\,t\right)\right\}$, or $\left[\left[u\left(x\,t\right)+t\,\mathrm{color}=\mathrm{blue}\right]\,w\left(x\,t\right)\right]$.

plot

The plot method supports the following calling sequences.

plot(v1=numeric, options)

plot(depspec, v1=numeric, options)

plot(v1=numeric, v2=range, options)

plot(depspec, v1=numeric, v2=range, options)

In the above calling sequences, the values/ranges specified by v1=numeric and v2=range, where v1,v2 are the names of the independent variables of the problem, specify the domain of the plot.

If only the value of v1 is specified (as is the case in the first two calling sequences), the functions are plotted with the value of that independent variable fixed as specified, and the domain of the plot ranging between the boundaries for the other independent variable. For example, specification of $t=1$ for a problem in $u\left(x\,t\right)$ defined for $x$ in the range $\left[0\,1\right]$ plots $u\left(x\,1\right)$ for $x=0..1$.

Note: If the spatial value is fixed, a time range must also be specified as in the third and fourth calling sequences.

If both a value and range are specified, v1 is fixed at the specified value and the plot is displayed over the range described by v2. The range for v2 must be contained in the boundaries for v2 (that is, the domain of the problem).

The options input specifies plot options to apply to all displayed plots that do not have a conflicting plot-specific option in the depspec.

Unlike plot, the plot method defaults to using the discrete points used on the computational domain, as controlled by timestep and/or spacestep. This can be changed by using the numpoints=npt plot option, in which case npt equally spaced points are computed by interpolation to obtain the plot.

animate

The animate method supports the following calling sequences.

animate(v1=numeric, range, list, options)

animate(depspec, v1=numeric, range, list, options)

animate(v1=numeric, range, list, v2=range, options)

animate(depspec, v1=numeric, range, list, v2=range, options)

In the above calling sequences, the values/ranges specified for v1 give the animation, and the range specified by v2=range specifies the domain of each frame of the animation where v1,v2 are the names of the independent variables of the problem.

The specification of the values for v1 determines which values of v1 correspond to each frame of the animation. If specified as a range, the default of 16 frames is used to generate the animation, which can be changed using the option frames=integer. If specified as a list, the v1 values in the list are sorted, and correspond to the v1 values for each frame of the animation. Note: The v1=numeric specification, where v1 is the time variable, corresponds to the specification $t=t_0..\mathrm{numeric}$ for $t_0<\mathrm{numeric}$ and $t=\mathrm{numeric}..t_0$ for $\mathrm{numeric}<t_0$, where $t_0$ is the initial time for the problem.

The v2=range argument (if provided) gives the range of the plot. If not present, this defaults to the interval between the boundaries of the problem.

Note: If v1 corresponds to the spatial variable, then v2=range must be provided.

The options input specifies plot options to apply to all displayed animations that do not have a conflicting plot-specific option in the depspec.

In addition, the title option has an additional feature. Using a floating point printf format code, you can display the time or space value for the current frame of the animation. For example, in a time-based animation, specifying title="t=%f" displays the specified title with the current time in place of %f in the animation. For more information, see the section on floating point formats on the printf help page.

Unlike plots[animate], the animate method defaults to using the discrete points used on the computational domain, as controlled by timestep or spacestep. This can be changed by using the numpoints=npt plot option, in which case npt equally spaced points are computed by interpolation to obtain the plot for each frame.

plot3d

The plot3d method supports the following calling sequences.

plot3d(v1=numeric, range, options)

plot3d(depspec, v1=numeric, range, options)

plot3d(v1=numeric, range, v2=range, options)

plot3d(depspec, v1=numeric, range, v2=range, options)

In the above calling sequences, the values/ranges specified for v1 and v2 specify the rectangular area over which the 3-d plot is produced where v1,v2 are the names of the independent variables of the problem.

The specifications of v1,v2 define the domain for the plot. Like animate, the specification v1=numeric defines the interval from the initial point of the problem for the time variable.

The options input specifies plot3d options to apply to all displayed plots that do not have a conflicting plot-specific option in the depspec.

Unlike plot3d, the plot3d method supports explicit functions of the dependent and independent variables in the color specification for both the RGB and HUE types. For example, for a PDE problem in $u\left(x\&comma;t\right)$, the RGB values can be specified by $\mathrm{color}=\left[u\left(x\&comma;t\right)\&comma;x-{u\left(x\&comma;t\right)}^2\&comma;\frac{xt}{u\left(x\&comma;t\right)}\right]$.

value

The value method provides a mechanism to obtain the numerical values for the solution of the PDE or PDE system directly. The value method supports the following calling sequences.

value(options)

value(depspec, options)

value(v1=numeric, options)

value(depspec, v1=numeric, options)

value(v1=numeric, v2=range, options)

value(depspec, v1=numeric, v2=range, options)

In all cases, the value method returns a procedure or procedures that can be used to evaluate the solution at a point.

There are two distinct forms of output, bivariate functions (that is, the solution over the domain), or univariate functions (that is, solution profiles with one independent variable value fixed).

The first two calling sequences provide bivariate function output, so the returned function or functions are procedures that, given the values of both independent variables, return solution value(s). For these calling sequences, the solution is computed as each point is requested.

The remaining four calling sequences provide solution profiles with v1 fixed to the specified numeric function value. For these calling sequences, the solution is computed prior to output, and the returned procedure essentially provides an interpolation of that output. For the 5th and 6th calling sequences, v1 corresponds to the spatial variable. In this case, the returned procedures can only compute values for the time range specified by v2=range.

The only allowable option is used to specify the output type for the returned procedure(s), and like dsolve[numeric], it can be either output=procedurelist (the default) or output=listprocedure. The procedurelist output returns a single procedure that, if called with valid independent variable values, returns a list of equations for the independent variable values, and the dependent variable values as specified by the depspec if provided. The listprocedure output returns a list of equations, where the right-hand sides of these equations correspond to procedures that can be called to compute the numerical value indicated by the left-hand sides of the equations.

settings

The settings command allows querying or setting of parameters for the computation of the solution. It supports the following calling sequences.

settings()

settings(options)

The first calling sequence is query-only, and reports settings that can be changed for the computation.

The second calling sequence changes settings based on the specified options. Valid options include timestep, and spacestep, which can be set in the same way as for the original call to pdsolve/numeric, with the exception that changing the spacestep value after creation of the solution module does not change the timestep. Adaptive modules have additional options for the settings method, see pdsolve,numeric,errorcontrol.

$\mathrm{PDE}≔\mathrm{diff}\left(u\left(x\&comma;t\right)\&comma;t\right)=-\mathrm{diff}\left(u\left(x\&comma;t\right)\&comma;x\right)$

$\mathrm{PDE}≔\frac{\partial }{\partial t}u\left(x\&comma;t\right)=-\frac{\partial }{\partial x}u\left(x\&comma;t\right)$

$\mathrm{IBC}≔\left\{u\left(0\&comma;t\right)=-\sin \left(2\mathrm{\pi }t\right)\&comma;u\left(x\&comma;0\right)=\sin \left(2\mathrm{\pi }x\right)\right\}$

$\mathrm{IBC}≔\left\{u\left(0\&comma;t\right)=-\sin \left(2\mathrm{\pi }t\right)\&comma;u\left(x\&comma;0\right)=\sin \left(2\mathrm{\pi }x\right)\right\}$

$\mathrm{pds}≔\mathrm{pdsolve}\left(\mathrm{PDE}\&comma;\mathrm{IBC}\&comma;\mathrm{numeric}\&comma;\mathrm{time}=t\&comma;\mathrm{range}=0..1\right)$

$\mathrm{pds}≔\mathbf{module}\left(\right)\mathbf{export}\mathrm{plot}\&comma;\mathrm{plot3d}\&comma;\mathrm{animate}\&comma;\mathrm{value}\&comma;\mathrm{settings}\&semi;...\mathbf{end\; module}$

$\mathrm{p1}≔\mathrm{pds}:-\mathrm{plot}\left(t=0\&comma;\mathrm{numpoints}=50\right)\&colon;$$\mathrm{p2}≔\mathrm{pds}:-\mathrm{plot}\left(t=\frac{1}{8}\&comma;\mathrm{numpoints}=50\&comma;\mathrm{color}=\mathrm{blue}\right)\&colon;$$\mathrm{p3}≔\mathrm{pds}:-\mathrm{plot}\left(t=\frac{1}{4}\&comma;\mathrm{numpoints}=50\&comma;\mathrm{color}=\mathrm{green}\right)\&colon;$$\mathrm{plots}\left[\mathrm{display}\right]\left(\left\{\mathrm{p1}\&comma;\mathrm{p2}\&comma;\mathrm{p3}\right\}\right)$

$\mathrm{esol}≔\sin \left(2\mathrm{\pi }\left(x-t\right)\right)$

$\mathrm{esol}≔\sin \left(2\mathrm{\pi }\left(x-t\right)\right)$

$\mathrm{p2}≔\mathrm{pds}:-\mathrm{plot}\left(u-\mathrm{esol}\&comma;t=\frac{1}{8}\&comma;\mathrm{numpoints}=50\&comma;\mathrm{color}=\mathrm{blue}\right)\&colon;$$\mathrm{p3}≔\mathrm{pds}:-\mathrm{plot}\left(u-\mathrm{esol}\&comma;t=\frac{3}{8}\&comma;\mathrm{numpoints}=50\&comma;\mathrm{color}=\mathrm{green}\right)\&colon;$$\mathrm{plots}\left[\mathrm{display}\right]\left(\left\{\mathrm{p2}\&comma;\mathrm{p3}\right\}\right)$

$\mathrm{PDE}≔\mathrm{diff}\left(u\left(x\&comma;t\right)\&comma;t\right)=\frac{1}{10}\mathrm{diff}\left(u\left(x\&comma;t\right)\&comma;x\&comma;x\right)$

$\mathrm{PDE}≔\frac{\partial }{\partial t}u\left(x\&comma;t\right)=\frac{\frac{{\partial }^2}{\partial x^2}u\left(x\&comma;t\right)}{10}$

$\mathrm{IBC}≔\left\{u\left(0\&comma;t\right)=0\&comma;u\left(x\&comma;0\right)=1\&comma;\mathrm{D}\left[1\right]\left(u\right)\left(1\&comma;t\right)=0\right\}$

$\mathrm{IBC}≔\left\{u\left(0\&comma;t\right)=0\&comma;u\left(x\&comma;0\right)=1\&comma;{\mathrm{D}}_1\left(u\right)\left(1\&comma;t\right)=0\right\}$

$\mathrm{pds}≔\mathrm{pdsolve}\left(\mathrm{PDE}\&comma;\mathrm{IBC}\&comma;\mathrm{numeric}\right)$

$\mathrm{pds}≔\mathbf{module}\left(\right)\mathbf{export}\mathrm{plot}\&comma;\mathrm{plot3d}\&comma;\mathrm{animate}\&comma;\mathrm{value}\&comma;\mathrm{settings}\&semi;...\mathbf{end\; module}$

$\mathrm{p1}≔\mathrm{pds}:-\mathrm{plot}\left(t=0\right)\&colon;$$\mathrm{p2}≔\mathrm{pds}:-\mathrm{plot}\left(t=\frac{1}{10}\right)\&colon;$$\mathrm{p3}≔\mathrm{pds}:-\mathrm{plot}\left(t=\frac{1}{2}\right)\&colon;$$\mathrm{p4}≔\mathrm{pds}:-\mathrm{plot}\left(t=1\right)\&colon;$$\mathrm{p5}≔\mathrm{pds}:-\mathrm{plot}\left(t=2\right)\&colon;$$\mathrm{plots}\left[\mathrm{display}\right]\left(\left\{\mathrm{p1}\&comma;\mathrm{p2}\&comma;\mathrm{p3}\&comma;\mathrm{p4}\&comma;\mathrm{p5}\right\}\&comma;\mathrm{title}=\mathrm{`Heat\; profile\; at\; t=0,0.1,0.5,1,2`}\right)$

$\mathrm{pds}:-\mathrm{value}\left(t=1\&comma;\mathrm{output}=\mathrm{listprocedure}\right)$

$\left[x=\mathbf{proc}\left(x\right)...\mathbf{end\; proc}\&comma;t=1.\&comma;u\left(x\&comma;t\right)=\mathbf{proc}\left(x\right)...\mathbf{end\; proc}\right]$

$\mathrm{uval}≔\mathrm{rhs}\left(\mathrm{op}\left(3\&comma;\right)\right)$

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

$\mathrm{fsolve}\left(\mathrm{uval}\left(x\right)=\frac{1}{2}\&comma;x=0..1\right)$

$0.2978753742$

$\mathrm{pds}:-\mathrm{plot3d}\left(t=0..1\&comma;x=0..1\&comma;\mathrm{axes}=\mathrm{boxed}\&comma;\mathrm{orientation}=\left[-120\&comma;40\right]\&comma;\mathrm{color}=\left[0\&comma;0\&comma;u\right]\right)$

$\mathrm{PDE}≔\mathrm{diff}\left(u\left(x\&comma;t\right)\&comma;t\&comma;t\right)=-\frac{1}{10}\mathrm{diff}\left(u\left(x\&comma;t\right)\&comma;x\&comma;x\&comma;x\&comma;x\right)$

$\mathrm{PDE}≔\frac{{\partial }^2}{\partial t^2}u\left(x\&comma;t\right)=-\frac{\frac{{\partial }^4}{\partial x^4}u\left(x\&comma;t\right)}{10}$

$\mathrm{IBC}≔\left\{u\left(1\&comma;t\right)=0\&comma;u\left(x\&comma;0\right)={\sin \left(\mathrm{\pi }x\right)}^2\&comma;u\left(-1\&comma;t\right)=0\&comma;\mathrm{D}\left[1\right]\left(u\right)\left(1\&comma;t\right)=0\&comma;\mathrm{D}\left[1\right]\left(u\right)\left(-1\&comma;t\right)=0\&comma;\mathrm{D}\left[2\right]\left(u\right)\left(x\&comma;0\right)=0\right\}$

$\mathrm{IBC}≔\left\{u\left(\mathrm{-1}\&comma;t\right)=0\&comma;u\left(1\&comma;t\right)=0\&comma;u\left(x\&comma;0\right)={\sin \left(\mathrm{\pi }x\right)}^2\&comma;{\mathrm{D}}_1\left(u\right)\left(\mathrm{-1}\&comma;t\right)=0\&comma;{\mathrm{D}}_1\left(u\right)\left(1\&comma;t\right)=0\&comma;{\mathrm{D}}_2\left(u\right)\left(x\&comma;0\right)=0\right\}$

$\mathrm{pds}≔\mathrm{pdsolve}\left(\mathrm{PDE}\&comma;\mathrm{IBC}\&comma;\mathrm{numeric}\&comma;\mathrm{spacestep}=\frac{1}{40}\right)$

$\mathrm{pds}≔\mathbf{module}\left(\right)\mathbf{export}\mathrm{plot}\&comma;\mathrm{plot3d}\&comma;\mathrm{animate}\&comma;\mathrm{value}\&comma;\mathrm{settings}\&semi;...\mathbf{end\; module}$

$\mathrm{pds}:-\mathrm{animate}\left(t=2\&comma;\mathrm{frames}=50\&comma;\mathrm{title}=''time\; =\; \%f''\right)$

$\mathrm{PDE}≔\mathrm{diff}\left(u\left(x\&comma;t\right)\&comma;t\right)=-\mathrm{diff}\left(u\left(x\&comma;t\right)\&comma;x\right)u\left(x\&comma;t\right)$

$\mathrm{PDE}≔\frac{\partial }{\partial t}u\left(x\&comma;t\right)=-\left(\frac{\partial }{\partial x}u\left(x\&comma;t\right)\right)u\left(x\&comma;t\right)$

$\mathrm{IBC}≔\left\{u\left(0\&comma;t\right)=\frac{\exp \left(-4\right)}{2}\&comma;u\left(x\&comma;0\right)=\frac{\exp \left(-16{\left(x-\frac{1}{2}\right)}^2\right)}{2}\right\}$

$\mathrm{IBC}≔\left\{u\left(0\&comma;t\right)=\frac{{\&ExponentialE;}^{\mathrm{-4}}}{2}\&comma;u\left(x\&comma;0\right)=\frac{{\&ExponentialE;}^{-16{\left(x-\frac{1}{2}\right)}^2}}{2}\right\}$

$\mathrm{pds}≔\mathrm{pdsolve}\left(\mathrm{PDE}\&comma;\mathrm{IBC}\&comma;\mathrm{numeric}\&comma;\mathrm{time}=t\&comma;\mathrm{range}=0..2\&comma;\mathrm{spacestep}=\frac{1}{256}\&comma;\mathrm{timestep}=\frac{1}{60}\right)$

$\mathrm{pds}≔\mathbf{module}\left(\right)\mathbf{export}\mathrm{plot}\&comma;\mathrm{plot3d}\&comma;\mathrm{animate}\&comma;\mathrm{value}\&comma;\mathrm{settings}\&semi;...\mathbf{end\; module}$

$\mathrm{pds}:-\mathrm{animate}\left(t=0.60\&comma;\mathrm{frames}=31\right)$

$\mathrm{PDE}≔\left\{\mathrm{diff}\left(u\left(x\&comma;t\right)\&comma;t\right)=-\frac{1}{20}\mathrm{diff}\left(v\left(x\&comma;t\right)\&comma;x\&comma;x\right)\&comma;\mathrm{diff}\left(v\left(x\&comma;t\right)\&comma;t\right)=\mathrm{diff}\left(u\left(x\&comma;t\right)\&comma;x\&comma;x\right)\right\}\&colon;$

$\mathrm{IBC}≔\left\{u\left(0\&comma;t\right)=0\&comma;u\left(1\&comma;t\right)=0\&comma;u\left(x\&comma;0\right)=\sin \left(\mathrm{\pi }x\right)\&comma;v\left(0\&comma;t\right)=1\&comma;v\left(1\&comma;t\right)=0\&comma;v\left(x\&comma;0\right)=1-x\right\}\&colon;$

$\mathrm{pds}≔\mathrm{pdsolve}\left(\mathrm{PDE}\&comma;\mathrm{IBC}\&comma;\mathrm{numeric}\&comma;\mathrm{spacestep}=\frac{1}{50}\right)$

$\mathrm{pds}≔\mathbf{module}\left(\right)\mathbf{export}\mathrm{plot}\&comma;\mathrm{plot3d}\&comma;\mathrm{animate}\&comma;\mathrm{value}\&comma;\mathrm{settings}\&semi;...\mathbf{end\; module}$

$\mathrm{plots}\left[\mathrm{display}\right]\left(\left[\mathrm{seq}\left(\mathrm{pds}:-\mathrm{plot}\left(v\&comma;t=\frac{i}{10}\right)\&comma;i=0..5\right)\right]\right)$

Boyce, W.E. and DiPrima, R.C. Elementary Differential Equations and Boundary Value Problems. New York: John Wiley & Sons, 1997.

Carrier, G.F. and Pearson, C.E. Partial Differential Equations: Theory and Technique. 2nd ed. Academic Press, 1988.

Strikwerda, J. C. Finite Difference Schemes and Partial Differential Equations. Wadsworth and Brooks/Cole, 1989.