We propose an algorithm based on variational quantum imaginary time evolution for solving the Feynman-Kac partial differential equation resulting from a multidimensional system of stochastic differential equations. We utilize the correspondence between the Feynman-Kac partial differential equation (PDE) and the Wick-rotated Schr{o}dinger equation for this purpose. The results for a $(2+1)$ dimensional Feynman-Kac system, obtained through the variational quantum algorithm are then compared against classical ODE solvers and Monte Carlo simulation. We see a remarkable agreement between the classical methods and the quantum variational method for an illustrative example on six qubits. In the non-trivial case of PDEs which are preserving probability distributions -- rather than preserving the $ell_2$-norm -- we introduce a proxy norm which is efficient in keeping the solution approximately normalized throughout the evolution. The algorithmic complexity and costs associated to this methodology, in particular for the extraction of properties of the solution, are investigated. Future research topics in the areas of quantitative finance and other types of PDEs are also discussed.