It can solve for scalars - and since both expressions are linear in the matrices, the scalar and the matrix solution coincide.
syms x0 x1 x2 m k p
eqn1 = m*x1*p^2+2*k*(x1-x0)+2*k*(x1-x2)==0;
eqn2 = m*x2*p^2+2*k*(x2-x1)==0;
sol = solve([eqn1,eqn2],[x1,x2])
sol.x1
sol.x2
