For the approximate solution of fractional partial differential equations, the authors suggest to use a method based on a (generalized) Taylor expansion of the solution. As in the classical case one can then try to compute the coefficients of the expansion, and hence the exact solution, by suitable recurrence relations based on the differential equation. This path seems to be viable only in the case of a sufficiently simple equation.
In particular, the method requires that the solution can be expanded in a series of a special form, and in a typical application it is by no means clear whether such an expansion is possible. Moreover it seems (see, e.g., Example 5.2) that the computed approximate solution does not depend on any boundary conditions. Obviously the exact solution does depend on the boundary conditions, and this discrepancy raises strong concerns about the correctness of the approach.