clear
reset(symengine)

syms k1 k2 positive

M = [
    0 0 1 0 ; 
    0 0 0 1 ; 
    -(k1^2 + k2^2) k2^2 0 0 ;
    k2^2 -(k1^2+k2^2) 0 0 ]

[U, D] = eig(M)


% rewrite U, D so that w1 = k1 and w2 = sqrt(k1^2 + 2*k2^2)
% that means that k1 = w1 and k2 = sqrt((w2^2 - w1^2)/2)

syms w1 w2 positive

k2expr = sqrt((w2^2 - w1^2)/2);

U = simple(subs(U, {k1,k2}, {w1, k2expr}))
D = simple(subs(D, {k1,k2}, {w1, k2expr}))

syms t real

B = U * expm(D*t)

syms x1 dx1 x2 dx2 real

q0 = [ x1 x2 dx1 dx2  ]'

alpha = inv(U) * q0

q_general = simplify(expand(B*alpha), 250)

q_case1 = simple(subs(q_general, {x2,dx1,dx2}, {x1,0,0}))
q_case2 = simple(subs(q_general, {x2,dx1,dx2}, {-x1,0,0}))
q_case3 = simple(subs(q_general, {x2,dx1,dx2}, {0,0,0}))