I am reading this paper1 for understanding of modified equation and numerical viscosity. In the first section of this paper it has many calculations involving with finite differencing, derivative, and series expansion. I am trying to use Mathematica to aid the symbolic computation. Here are some notes.

The Series function is used for series expansion. And the Normal function converts the series data into a polynomial.

For single variable function we can use following to get an second order series expansion at $$x0$$:

This answer from Mathematica.StackExchange shows how to find the expansion of a multi-variable function $$u$$ at $$(x0, t0)$$:

Now we can substitute the expanded function into the differencing scheme. Let us take the second-order Lax-Wendroff scheme as an example (the equation 1.3 in the paper):

where $u_j^n$, $u_j^{n+1}$, $u_{j+1}^n$ and $u_{j-1}^n$ are substituted by

Be aware that Superscript function should be used explicitly to denote a superscript. Otherwise terms such as $$u_j^n$$ would be recognized as exponents if they were input in the Ctrl ^ way.

The final result(the equation 1.5 in the paper) is

This result is rendered by a tool called pdConv which is introduced by Mathematica Q&A Series.

1. Warming, R. F., and B. J. Hyett. “The modified equation approach to the stability and accuracy analysis of finite-difference methods.” Journal of computational physics 14.2 (1974): 159-179.