This paper describes a technique for computing optimal upper bound limit load solution, based on the linear matching method, for the bearing capacity of surface strip footings on multilayer soils. The method assumes a perfectly plastic soil model with a Drucker–Prager yield criterion and involves the sequential matching of an appropriate linear material to the yield function. In this paper, results are presented in the form of bearing capacity factors based on various layer properties and geometries. Solutions for both an associated and a non-associated flow rules are discussed. Comparisons with existing analytic and finite element solutions demonstrate the efficiency of the current method.