Repository: sympy/sympy

Issue:
system om 2nd order ode
I use sympy to solve this problem. MWE:
```python
from sympy import *
import numpy as np
import matplotlib.pyplot as plt
x1,x2=symbols('x1,x2',cls=Function)
t,k1,k2,m1,m2=symbols('t k1 k2 m1 m2',real=True)
k1=1
k2=1
m1=1
m2=1
eq1=Eq(m1*diff(x1(t),t,2)+k1*x1(t)-k2*(x2(t)-x1(t)),0)  
eq2=Eq(m2*diff(x2(t),t,2)+k2*(x2(t)-x1(t)),0)  
eqs=[eq1,eq2]
IC1={x1(0):0,x1(t).diff().subs(t,0):0}
IC2={x2(0):1,x2(t).diff().subs(t,0):0}
[x1sol,x2sol]=dsolve(eqs,[x1(t),x2(t)],ics={x1(0):0,x1(t).diff().subs(t,0):0,x2(0):1,x2(t).diff().subs(t,0):0})
x1sol=simplify(x1sol)
print(x1sol)
x2sol=simplify(x2sol)
print(x2sol)
```
However, the solution for x1 and x2 doesn't fit the initial conditions. Maybe a bug in sympy? I report this as an issue.

Discussion:

Patch:
diff --git a/sympy/solvers/ode/ode.py b/sympy/solvers/ode/ode.py
index ac24403d5f..1da6ac9cff 100644
--- a/sympy/solvers/ode/ode.py
+++ b/sympy/solvers/ode/ode.py
@@ -770,10 +770,9 @@ def solve_ics(sols, funcs, constants, ics):
                 x0 = funcarg.variables[0]
                 variables = (x,)*len(funcarg.variables)
                 matching_func = deriv.subs(x0, x)
-            if variables not in diff_variables:
-                for sol in sols:
-                    if sol.has(deriv.expr.func):
-                        diff_sols.append(Eq(sol.lhs.diff(*variables), sol.rhs.diff(*variables)))
+            for sol in sols:
+                if sol.has(deriv.expr.func):
+                    diff_sols.append(Eq(sol.lhs.diff(*variables), sol.rhs.diff(*variables)))
             diff_variables.add(variables)
             S = diff_sols
         else:


Test patch:
diff --git a/sympy/solvers/ode/tests/test_ode.py b/sympy/solvers/ode/tests/test_ode.py
index 354c1073c2..c53d38d38b 100644
--- a/sympy/solvers/ode/tests/test_ode.py
+++ b/sympy/solvers/ode/tests/test_ode.py
@@ -1045,3 +1045,17 @@ def test_issue_22523():
         5.33333333333333*s**4*sqrt(N - 1.0)/N**3) + C1*s*(1.0 -
         1.33333333333333*s**3*sqrt(N - 1.0)/N - 0.666666666666667*s**2*sqrt(N
         - 1.0)/N + 1.33333333333333*s**3*sqrt(N - 1.0)/N**2) + O(s**6))'''))
+
+
+def test_issue_22604():
+    x1, x2 = symbols('x1, x2', cls = Function)
+    t, k1, k2, m1, m2 = symbols('t k1 k2 m1 m2', real = True)
+    k1, k2, m1, m2 = 1, 1, 1, 1
+    eq1 = Eq(m1*diff(x1(t), t, 2) + k1*x1(t) - k2*(x2(t) - x1(t)), 0)
+    eq2 = Eq(m2*diff(x2(t), t, 2) + k2*(x2(t) - x1(t)), 0)
+    eqs = [eq1, eq2]
+    [x1sol, x2sol] = dsolve(eqs, [x1(t), x2(t)], ics = {x1(0):0, x1(t).diff().subs(t,0):0, \
+                                                        x2(0):1, x2(t).diff().subs(t,0):0})
+    assert x1sol == Eq(x1(t), sqrt(3 - sqrt(5))*(sqrt(10) + 5*sqrt(2))*cos(sqrt(2)*t*sqrt(3 - sqrt(5))/2)/20 + \
+                       (-5*sqrt(2) + sqrt(10))*sqrt(sqrt(5) + 3)*cos(sqrt(2)*t*sqrt(sqrt(5) + 3)/2)/20)
+    assert x2sol == Eq(x2(t), (sqrt(5) + 5)*cos(sqrt(2)*t*sqrt(3 - sqrt(5))/2)/10 + (5 - sqrt(5))*cos(sqrt(2)*t*sqrt(sqrt(5) + 3)/2)/10)