We describe a numerical method for solving 3-dimensional partial differential equations, which arise in mathematical finance and other applications. The goal of the paper is to introduce a technique based on Wiener-Hopf factorization with application of Laplace transform. We analyze the problem in terms of expectations of random processes. We construct an approximation scheme by using Carr randomization and constructing a Markov chain, and reduce the original problem to a sequence of 1-dimensional integro-differential equations with suitable boundary conditions. The kernels of the equations are defined by Levy processes with constant variance. An analytic solution to each problem can be expressed in terms of Laplace-Carson transform of the corresponding characteristic functions of its supremum and infimum processes. We show that for a class of models it is possible to construct an efficient method for solving these equations which relies upon approximate formulae for the transform, and discuss modifications allowing to reduce the amount of computations.