The present study establishes an accurate and efficient algorithm based on Monte Carlo (MC) simulation for solving high dimensional linear systems of algebraic equations (LSAEs) and two-dimensional Fredholm integral equations of the second kind (FIESK). This new combined numerical-probabilistic algorithm is based on Jacobi over-relaxation method and MC simulation in conjunction with the iterative refinement technique to find the unique solution of the large sparse LSAEs. It has an excellent accuracy, low cost and simple structure. Theoretical results are established to justify the convergence of the algorithm. To confirm the accuracy and efficiency of the present work, the proposed algorithm is used for solving 1000 £ 1000 and 2000 £ 2000 LSAEs. Furthermore, the algorithm is coupled with Galerkin's method to illustrate the power and effectiveness of the proposed algorithm for solving twodimensional FIESK.