This study focuses on the topic of model order reduction (MOR) for coupled systems with inhomogeneous initial conditions and presents an MOR method by general orthogonal polynomials with Arnoldi algorithm. The main procedure is to use a series of expansion coefficients vectors in the space spanned by orthogonal polynomials that satisfy a recursive formula to generate a projection based on the multiorder Arnoldi algorithm. The resulting model not only match desired number of expansion coefficients but also has the same coupled structure as the original system. Moreover, the stability is preserved as well. The error bound between the outputs is well‐discussed. Finally, numerical results show that the authors’ method can deal well with those systems with inhomogeneous initial conditions in the views of accuracy and computational cost.