Fractional order integrodifferential equations cannot be directly solved like ordinary differential equations. Numerical methods for such equations have additional algorithmic complexities. We present a particularly simple recipe for solving such equations using a Galerkin scheme developed in prior work. In particular, matrices needed for that method have here been precisely evaluated in closed form using special functions, and a small Matlab program is provided for the same. For equations where the highest order of the derivative is fractional, differential algebraic equations arise; however, it is demonstrated that there is a simple regularization scheme that works for these systems, such that accurate solutions can be easily obtained using standard solvers for stiff differential equations. Finally, the role of nonzero initial conditions is discussed in the context of the present approximation method.