Two new methods for the solution of problems involving material variability are proposed. Medium properties are modeled as second order stochastic processes defined by their mean and convariance functions. Both methods make use of the Karhunen-Loeve expansion which is a mean-square convergent orthogonal expansion of a continuous process in terms of a countable set of uncorrelated random variables. The first of the proposed methods relies on implementing the Karhunen-Loeve expansion for the medium property in conjunction with a Neumann expansion of the inverse operator. This results in an explicit expression for the response process as a multivariate polynomial functional of a set of uncorrelated random variables. The second method treats the solution process as an element in the Hilbert space of random functions, in which a sequence of projection operators is identified as the Polynomial Chaoses of consecutive order. The solution process is then determined by its projection onto the spaces spanned by these polynomials. These concepts can be construed as extensions of the deterministic finite element methods to the space of random functions. Both of the proposed methods are exemplified by three problems from the field of engineering mechanics. The corresponding results are found in agreement with those obtained by a Monte-Carlo simulation solution of the problems. In addition to the two methods mentioned above, a new formulation is presented for a class of problems involving deterministic media subjected to random external excitations. The formulation involves a combination of the boundary element method with the Karhunen-Loeve expansion for the exciting process. Namely, the boundary element method is used as a discretization tool to restate the problem as a set of discrete equations. Further, the Karhunen-Loeve expansion is utilized to represent the random processes in a manner conducing their optimal discretization.