This thesis deals with developing Galerkin based solution strategies for several important classes of differential equations involving derivatives and integrals of various fractional orders. Fractional order calculus finds use in several areas of science and engineering. The use of fractional derivatives may arise purely from the mathematical viewpoint, as in controller design, or it may arise from the underlying physics of the material, as in the damping behavior of viscoelastic materials. The physical origins of the fractional damping motivated us to study viscoelastic behavior of disordered materials at three levels. At the first level, we review two first principles models of rubber viscoelasticity. This leads us to study, at the next two levels, two simple disordered systems. The study of these two simplified systems prompted us towards an infinite dimensional system which is mathematically equivalent to a fractional order derivative or integral. This infinite dimensional system forms the starting point for our Galerkin projection based approximation scheme. In a simplified study of disordered viscoelastic materials, we show that the networks of springs and dash-pots can lead to fractional power law relaxation if the damping coefficients of the dash-pots follow a certain type of random distribution. Similar results are obtained when we consider a more simplified model, which involves a random system coefficient matrix. Fractional order derivatives and integrals are infinite dimensional operators and non-local in time: the history of the state variable is needed to evaluate such operators. This non-local nature leads to expensive long-time computations (O(t2) computations for solution up to time t). A finite dimensional approximation of the fractional order derivative can alleviate this problem. We present one such approximation using a Galerkin projection. The original infinite dimensional system is replaced with an equivalent infinite dimensional system involving a partial differential equation (PDE). The Galerkin projection reduces the PDE to a finite system of ODEs. These ODEs can be solved cheaply (O(t) computations). The shape functions used for the Galerkin projection are important, and given attention. Calculations with both global shape functions as well as finite elements are presented. The discretization strategy is improved in a few steps until, finally, very good performance is obtained over a user-specifiable frequency range (not including zero). In particular, numerical examples are presented showing good performance for frequencies varying over more than 7 orders of magnitude. For any discretization held fixed, however, errors will be significant at sufficiently low or high frequencies. We discuss why such asymptotics may not significantly impact the engineering utility of the method. Following this, we identify eight important classes of fractional differential equations (FDEs) and fractional integrodifferential equations (FIEs), and develop separate Galerkin based solution strategies for each of them. Distinction between these classes arises from the fact that both Riemann-Liouville as well as Caputo type derivatives used in this work do not, in general, follow either the law of exponents or the commutative property. Criteria used to identify these classes include; the initial conditions used, order of the highest derivative, integer or fractional order highest derivative, single or multiterm fractional derivatives and integrals. A key feature of our approximation scheme is the development of differential algebraic equations (DAEs) when the highest order derivative is fractional or the equation involves fractional integrals only. To demonstrate the effectiveness of our approximation scheme, we compare the numerical results with analytical solutions, when available, or with suitably developed series solutions. Our approximation scheme matches analytical/series solutions very well for all classes considered.