A difference equation simulator for a linear physical system specified in terms of either a differential equation or a rational function ofscan be derived by applying the bilinear transform. It is shown how the method may be extended to enable simulation of nonlinear systems. The same formal procedure can be used with higher‐order (e.g., Runge–Kutta) integrators to give considerably enhanced accuracy for the same step size. Three simple examples are given illustrating the method.