  A new method that enables easy and convenient discretization of partial
differential equations with derivatives of arbitrary real order (so-called
fractional derivatives) and delays is presented and illustrated on numerical
solution of various types of fractional diffusion equation. The suggested
method is the development of Podlubny's matrix approach (Fractional Calculus
and Applied Analysis, vol. 3, no. 4, 2000, 359--386). Four examples of
numerical solution of fractional diffusion equation with various combinations
of time/space fractional derivatives (integer/integer, fractional/integer,
integer/fractional, and fractional/fractional) with respect to time and to the
spatial variable are provided in order to illustrate how simple and general is
the suggested approach. The fifth example illustrates that the method can be
equally simply used for fractional differential equations with delays. A set of
MATLAB routines for the implementation of the method as well as sample code
used to solve the examples have been developed.
