A method is presented which can be used to obtain analytical solutions for boundary value problems (BVPs) using the matrix exponential and Maple. Systems of second order, linear differential equations are expressed as two or more first order equations in matrix form, and their solutions are obtained using the matrix exponential, matrix integration, and matrix inverse methods using Maple. The solution process is illustrated for single and multiple domains with different types of boundary conditions and constraints when necessary due to the boundary conditions. The method is easier to use and could be extended to include partial differential equations (PDEs).