The forward time forward/backward space method is an explicit single stage method that can be used to find solutions to first order time/space PDEs.

The method is $\mathrm{O}\left(h\,k\right)$ accurate.

First order PDEs that describe right-traveling waves should use the backward method (specified as ForwardTime1Space[backward]), and specify the boundary condition at the left boundary, while left-traveling waves require the forward method (specified as ForwardTime1Space[forward]), and the boundary condition at the right boundary. Numerical boundary conditions are not required.

As this is an explicit scheme, stability requires a restriction of the form $k<ah$ for some $a$ depending upon the problem.

The scheme used to compute the value at the mesh point $\left[1\&comma;i\right]$ is the result of applying differencing to the PDE about the point $t,x_i$ and can obtained by substitution of the discretizations

into the PDE for the backward method. The forward method simply replaces the discretization for $u_x$ by

The centered time forward/backward space method is an implicit single stage method that can be used to find solutions to PDEs containing the derivatives $u,u_x,u_t,u_{\mathrm{xt}}$.

The method is $\mathrm{O}\left(h\&comma;k^2\right)$ accurate.

PDEs that describe right-traveling waves should use the backward method (specified as CenteredTime1Space[backward]), and specify the boundary condition at the left boundary, while left-traveling waves require the forward method (specified as CenteredTime1Space[forward]), and the boundary condition at the right boundary. Numerical boundary conditions are not required.

This is an implicit scheme, and is unconditionally stable for many problems (though this may need to be checked).

The scheme used to compute the value at the mesh point $\left[1\&comma;i\right]$ is the result of applying differencing to the PDE about the point $t+\frac{k}{2},x_i$ and is obtained by substitution of the discretizations

into the PDE for the backward method. The forward method simply replaces the discretizations for $u_{\mathrm{xt}}$ and $u_x$ by

The backward time forward/backward space method is an implicit single stage method that can be used to find solutions to PDEs containing the derivatives $u,u_x,u_t,u_{\mathrm{xt}}$.

The method is $\mathrm{O}\left(h\&comma;k\right)$ accurate.

PDEs that describe right-traveling waves should use the backward method (specified as BackwardTime1Space[backward]), and specify the boundary condition at the left boundary, while left-traveling waves require the forward method (specified as BackwardTime1Space[forward]), and the boundary condition at the right boundary. Numerical boundary conditions are not required.

This is an implicit scheme, and is unconditionally stable for many problems (though this may need to be checked).

The scheme used to compute the value at the mesh point $\left[1\&comma;i\right]$ is the result of applying differencing to the PDE about the point $t+k,x_i$ and is obtained by substitution of the discretizations

into the PDE for the backward method. The forward method simply replaces the discretizations for $u_{\mathrm{xt}}$ and  $u_x$ by

The forward time centered space method is an explicit single stage method that can be used to find solutions to PDEs that are first order in time, and arbitrary order in space, with no mixed partial derivatives.

The method is $\mathrm{O}\left(h^2\&comma;k\right)$ accurate.

PDEs that are odd order in space require that numerical boundary conditions be specified so that each boundary has the same total number of conditions, and that the total number of boundary conditions is one greater than the differential order in the space variable. PDEs that are second order in space require no numerical boundary conditions, but still require the same number of conditions on each boundary.

As this is an explicit scheme, stability imposes a restriction on the time step of the problem relative to the space step. In general, this scheme is unstable for PDEs where $\mathrm{xord}$ is odd. For PDEs that are even order in space, stability requires a restriction of the form $k<a{h}^{\mathrm{xord}}$ for some $a$ depending upon the problem.

The scheme used to compute the value at the mesh point $\left[1\&comma;i\right]$ is the result of applying central differencing to the PDE about the point $t,x_i$ using $1+2⌈\frac{\mathrm{xord}}{2}⌉$ points (that is, the mesh uses $⌈\frac{\mathrm{xord}}{2}⌉$ points on either side of the center).

The centered time centered space method is an implicit single stage method that can be used to find solutions to PDEs that are first order in time, and arbitrary order in space, with mixed partial derivatives.

The method is $\mathrm{O}\left(h^2\&comma;k^2\right)$ accurate.

PDEs that are odd order in space require that numerical boundary conditions be specified so that each boundary has the same total number of conditions, and that the total number of boundary conditions is one greater than the differential order in the space variable. PDEs that are second order in space require no numerical boundary conditions, but still require the same number of conditions on each boundary.

This is an implicit scheme, and is unconditionally stable for many problems (though this may need to be checked).

The scheme used to compute the value at the mesh point $\left[1\&comma;i\right]$ is the result of applying central differencing to the PDE about the point $t+\frac{k}{2},x_i$ using $1+2⌈\frac{\mathrm{xord}}{2}⌉$ points (that is, the mesh uses $⌈\frac{\mathrm{xord}}{2}⌉$ points on either side of the center).

The backward time centered space method is an implicit method that can be used to find solutions to PDEs that are first order in time, and arbitrary order in space, with mixed partial derivatives.

The method is $\mathrm{O}\left(h^2\&comma;k\right)$ accurate.

PDEs that are odd order in space require that numerical boundary conditions be specified so that each boundary has the same total number of conditions, and that the total number of boundary conditions is one greater than the differential order in the space variable. PDEs that are second order in space require no numerical boundary conditions, but still require the same number of conditions on each boundary.

This is an implicit scheme, and is unconditionally stable for many problems (though this may need to be checked).

The scheme used to compute the value at the mesh point $\left[1\&comma;i\right]$ is the result of applying central differencing to the PDE about the point $t+k,x_i$ using $1+2⌈\frac{\mathrm{xord}}{2}⌉$ points (that is, the mesh uses $⌈\frac{\mathrm{xord}}{2}⌉$ points on either side of the center). For PDEs having no mixed partial derivatives, the only point in the discretization that is at the current time step is $u_{0,i}$, all other points are of the form $u\&lsqb;1,*\&rsqb;$.

The box method is an implicit single stage method that can be used to find solutions to PDEs that are first order in time, and arbitrary order in space, with mixed partial derivatives.

The method is $\mathrm{O}\left(h^2\&comma;k^2\right)$ accurate.

The method uses a spatial mesh with an even number of points, so necessarily one boundary must have one condition more than the other. Since the method is symmetric (would be no different if one boundary had the extra condition versus the other) the method is automatically adjusted to accommodate the boundary conditions.

PDEs that are even order in space require one additional numerical boundary condition, while PDEs that are odd order in space require no numerical boundary conditions.

This is an implicit scheme, and is unconditionally stable for many problems (though this may need to be checked).

The scheme used to compute the value at the mesh point $\left[1\&comma;i\right]$ is different depending if the extra boundary condition (or numerical boundary condition for even order problems) is given at the left or right boundaries. If the extra condition is on the left, the scheme can be obtained through application of central differencing of the PDE about the point $t+\frac{k}{2},x_i-\frac{h}{2}$ using $1+⌊\frac{\mathrm{xord}}{2}⌋$ points to the left of $i$, and $⌊\frac{\mathrm{xord}}{2}⌋$ points to the right. If the extra condition is on the right, the scheme can be obtained through application of central differencing of the PDE about the point $t+\frac{k}{2},x_i+\frac{h}{2}$ using $⌊\frac{\mathrm{xord}}{2}⌋$ points to the left of $i$, and $1+⌊\frac{\mathrm{xord}}{2}⌋$ points to the right.

The Lax-Friedrichs method is an explicit single stage method that can be used to find solutions to PDEs that are first order in time, and odd order in space, with no mixed partial derivatives.

The method is $\mathrm{O}\left(h^2\&comma;k\&comma;\frac{h^{2⌈\frac{\mathrm{xord}}{2}⌉}}{k}\right)$ accurate. Note: The method can be applied to PDEs with odd order in space, but the stability requirement (see below) forces the error to be $\mathrm{O}\left(1\right)$ due to the $\mathrm{O}\left(\frac{h^{2⌈\frac{\mathrm{xord}}{2}⌉}}{k}\right)$ term.

PDEs that are odd order in space require that numerical boundary conditions be specified so that each boundary has the same total number of conditions, and that the total number of boundary conditions is one greater than the differential order in the space variable. PDEs that are second order in space require no numerical boundary conditions, but do require the same number of conditions on each boundary.

As this is an explicit scheme, stability imposes a restriction on the time step of the problem relative to the space step. For PDEs that are odd order in space, stability requires a restriction of the form $k<a{h}^{\mathrm{xord}}$ for some $a$ depending upon the problem. When this condition is combined with the accuracy of the method, the error is discovered to be $\mathrm{O}\left(h^2,a\&apos;h^{\mathrm{xord}},\frac{h}{a\&apos;}\right)$, or simply $\mathrm{O}\left(h\right)$ when $k=a\&apos;h^{\mathrm{xord}},a\&apos;<a$.

The scheme used to compute the value at the mesh point $\left[1\&comma;i\right]$ is identical to that of the Euler scheme with the exception that the value $u_{0,i}$ is replaced by the interpolation of the value at $u_{0,i}$ using all other points in the Euler stencil $u_{0,j}$ where $j\ne i$.

The Lax-Wendroff method is an explicit single stage method that can be used to find solutions to linear PDEs that are first order in time and space.

The method is $\mathrm{O}\left(h^2\&comma;kh\&comma;k^2\right)$ accurate.

The Lax-Wendroff method always requires a numerical boundary condition so that one boundary condition is specified for each boundary.

As this is an explicit scheme, stability requires a restriction of the form $k<ah$ for some $a$ depending upon the problem.

The derivation of this scheme is complex, and involves computing difference schemes of derivatives of the input PDE, (which results in the linearity restriction), but is available in many standard textbooks.

The Leapfrog method is an explicit two stage method that can be used to find solutions to PDEs that are first order in time, and arbitrary order in space, with no mixed partial derivatives.

The method is $\mathrm{O}\left(h^2\&comma;k^2\right)$ accurate, and as it is a two stage method, a startup method must be used.

PDEs that are odd order in space require that numerical boundary conditions be specified so that each boundary has the same total number of conditions, and that the total number of boundary conditions is one greater than the differential order in the space variable. PDEs that are even order in space require no numerical boundary conditions, but still require the same number of conditions on each boundary.

As this is an explicit scheme, stability imposes a restriction on the time step of the problem relative to the space step. In general, this scheme is unstable for PDEs where $\mathrm{xord}$ is even. For PDEs that are odd order in space, stability requires a restriction of the form $k<a{h}^{\mathrm{xord}}$ for some $a$ depending upon the problem.

The scheme used to compute the value at the mesh point $\left[1\&comma;i\right]$ uses the formula

for the time derivative, and uses central differencing over $1+2⌈\frac{\mathrm{xord}}{2}⌉$ points to compute all spatial derivatives.

The DuFort-Frankel method is an explicit two stage method that can be used to find solutions to linear PDEs that are first order in time, and even order in space, with no mixed partial derivatives. In some cases it can be used to find solutions of nonlinear PDEs.

The method is $\mathrm{O}\left(h^2\&comma;k^2\&comma;\frac{k^2}{h^{2⌈\frac{\mathrm{xord}}{2}⌉}}\right)$ accurate, and as it is a two stage method, a startup method must be used. Note: The method can be applied to PDEs with odd order in space, but the stability requirement (see below) forces the error to be $\mathrm{O}\left(1\right)$ due to the $\mathrm{O}\left(\frac{k^2}{h^{2⌈\frac{\mathrm{xord}}{2}⌉}}\right)$ term. Note, however, that if no even order derivatives appear in the PDE at all, then the resulting scheme is identical to the scheme generated by the Leapfrog method for the same PDE, so in that case those comments apply.

Since this method is restricted to PDEs that are even order in space, no numerical boundary conditions are required, but the number of conditions on each boundary must be the same.

As this is an explicit scheme, stability imposes a restriction on the time step of the problem relative to the space step. The restriction is of the form $k<a{h}^{\mathrm{xord}}$ for some $a$ depending upon the problem. When this condition is combined with the accuracy of the method, the error is discovered to be $\mathrm{O}\left(h^2,a{\&apos;}^2{h}^{\mathrm{xord}}\right)$, or simply $\mathrm{O}\left(h^2\right)$ when $k=a\&apos;h^{\mathrm{xord}},a\&apos;<a$.

The scheme used to compute the value at the mesh point $\left[1\&comma;i\right]$ is identical to that of the Leapfrog scheme with the exception that the value $u_{0,i}$ is replaced by the average of the values $u_{1,i}$ and $u_{\mathrm{-1},i}$ for all occurrences in the scheme.

Construction of a new method or scheme for use by pdsolve is somewhat involved, and requires some basic Maple programming knowledge. In essence, you must construct a routine that derives the desired method for an arbitrary point $\left[1\&comma;\mathrm{\xi }\right]$, for an input PDE, and returns that information, and information on the method itself in the entries of an input table.

Construction of new methods has the following limitations.

- The PDE must be first order in time.

- The scheme cannot use more than two additional stages.

- The time and space step sizes are constant.

- The scheme cannot use discrete values of the space variable that are outside the stencil of the method.

The following is the code for the CenteredTime1Space method and is used as an example to describe what is needed.

Input Data

The input data is fairly well described by the comments present in the parameter list of the example, only a few items need be discussed.

The name of the procedure for the method must be of the form `pdsolve/numeric/findif/<method name>`, as pdsolve checks for the presence of the routine based on the argument of the 'method' option.

The input PDE is an expression (not an equation) that is equivalent to the PDE being numerically solved by pdsolve. If the input to pdsolve is an equation, it is converted to an expression. It is converted to diff notation (not operator notation), so if ${\mathrm{D}}_1\left(u\right)\left(x\&comma;t\right)$ is present in the input PDE, this routine sees $\frac{\partial }{\partial x}u\left(x\&comma;t\right)$ instead.

The $'v'$ input is the base name that should be used to construct the indexed names for the discretization, that is, the value at the mesh point $\left[1\&comma;\mathrm{ix}\right]$ should be denoted as $v_{1,\mathrm{ix}}$.

The returned expression is to be used to compute the discrete solution value at $v_{1,\mathrm{ix}}$, where $'\mathrm{ix}'$ is the input parameter.

The INFO table has entries that can be used to more easily determine if the method being derived is valid for the input PDE and is also used to output certain information about the discretization. In addition, the $''directional''$ flag can be set as output.

The INFO table has the following entries that can be used by this routine.

"timeord"

The INFO["timeord"] entry holds an integer giving the differential order of the PDE with respect to the time variable.

"spaceord"

The INFO["spaceord"] entry holds an integer giving the differential order of the PDE with respect to the space variable.

"mixed"

The INFO["mixed"] entry holds boolean value indicating if mixed partial derivatives are present in the PDE.

"directional"

The INFO["directional"] flag should be set to 'true' by the routine for any implicit methods where the specific solving unknown is important. Otherwise, numeric may shift the method depending on the location of the boundary conditions. If not set, it is assumed to be false.

For example, the Box method applied to a first order PDE is of the same form if the boundary condition is given on the left, as if the boundary condition is given on the right. As a result, for this scheme, INFO["directional"] should be false (the default).

Note: Explicit methods have INFO["directional"] set to true automatically, as there is only ever one discrete dependent variable value present in the scheme for a future time step.

The rest is fairly straightforward. The routine should generate a difference scheme for the PDE, and the only unknowns present in the scheme should be the value of the dependent variable at mesh points, $v_{\mathrm{timeidx},\mathrm{spaceidx}}$, the step sizes $h,k$, the time value at the 0 time step $t$, and the space value on the mesh $x_{\mathrm{spaceidx}}$.

Additional Notes

The value of $'t'$ used in the application of the output scheme is as at the $0$ time step, that is, if the discretization is to be used to compute the discrete solution values at $t=\mathrm{tv}+k$, the time value $t$ used in computing with the scheme would be $\mathrm{tv}$. This means that schemes that are derived as differences at half mesh points or backward schemes that use the time value at the $'\mathrm{future}'$ solution point must adjust the $t$ value in the output scheme. In the latter case this can be accomplished by setting $t$ to $t+k$.

The pdsolve command performs well for a direct discretization of the PDE, due to the optimization settings, so use of $\mathrm{expand}$ or $\mathrm{normal}$ on the discretized PDE within the method procedure can degrade performance. It is best to experiment with this when writing a custom method.