This paper presents an algorithm to solve linear systems expressed by a matrix stored in a tensor product format. The proposed solution is based on a LU decomposition of the matrix keeping the tensor product structure. It is shown that the complexity of the decomposition is negligible and the backward and forward substitutions are no more complex than two standard vector-matrices multiplications. Finally, applications of the proposed algorithm and the comparison with other similar techniques are discussed.