Let $G=(V,E)$ be an undirected graph, $L_G\in \mathbb{R}^{V \times V}$ be the associated Laplacian matrix, and $b \in \mathbb{R}^V$ be a vector. Solving the Laplacian system $L_G x = b$ has numerous applications in theoretical computer science, machine learning, and network analysis. Recently, the notion of the Laplacian operator $L_F:\mathbb{R}^V \to 2^{\mathbb{R}^V}$ for a submodular transformation $F:2^V \to \mathbb{R}_+^E$ was introduced, which can handle undirected graphs, directed graphs, hypergraphs, and joint distributions in a unified manner. In this study, we show that the submodular Laplacian system $L_F( x) \ni b$ can be solved in polynomial time. Furthermore, we also prove that even when the submodular Laplacian system has no solution, we can solve its regression form in polynomial time. Finally, we discuss potential applications of submodular Laplacian systems in machine learning and network analysis.