> restart; with(LinearAlgebra); assume(omega, real, omega > 0);
> G := 9;
> z := (xi^2+xi/(1+xi^2))/(1+xi^2);
`output redirected...`> print(); # input placeholder
> C := `<,>`(1-z, seq(sin((n-1)*Pi*z), n = 2 .. G));
`output redirected...`> print(); # input placeholder
> g := Transpose(C);
`output redirected...`> print(); # input placeholder
> A := Multiply(C, g);
`output redirected...`> print(); # input placeholder
> B := xi^2*A;
`output redirected...`> print(); # input placeholder
> for nn to G do for hh to G do A[nn, hh] := evalf(int(A[nn, hh], xi = 0 .. infinity)); B[nn, hh] := evalf(int(B[nn, hh], xi = 0 .. infinity)) end do end do;
> for nn to G do C(nn, 1) := evalf(int(C(nn, 1), xi = 0 .. infinity)) end do;
> f := Transpose(C);
> H := MatrixInverse(A);
`output redirected...`> print(); # input placeholder
> M := Multiply(H, B);
`output redirected...`> print(); # input placeholder
> N := Multiply(H, C);
> gh := `<|>`(N, Vector(1 .. G, 0), -M);
`output redirected...`> print(); # input placeholder
> fgh2 := `<|>`(h1, K1-1, (7*K2*(1/11))*f);
> cv := `<|>`(1, Vector[row](1 .. G+1, 0));
`output redirected...`> print(); # input placeholder
> capB := `<,>`(fgh2, cv, gh);
`output redirected...`> print(); # input placeholder
> Sb := CharacteristicMatrix(capB, s);
`output redirected...`> print(); # input placeholder
> carpol2 := LinearAlgebra[Determinant](Sb);
`output redirected...`> print(); # input placeholder
> eq2 := subs(s = I*omega, carpol2);
`output redirected...`> print(); # input placeholder
> eq21 := evalc(Re(eq2));
`output redirected...`> print(); # input placeholder
> eq22 := evalc(Im(eq2));
`output redirected...`> print(); # input placeholder
> ee2 := solve({eq21, eq22}, {K1, K2});
`output redirected...`> print(); # input placeholder
> hio := subs(ee2[1], ee2[2], fgh2);
`output redirected...`> print(); # input placeholder
> capA := `<,>`(hio, cv, gh);
> kv := `<,>`(1, seq(v[nn], nn = 2 .. G+2));
`output redirected...`> print(); # input placeholder
> Iden := IdentityMatrix(G+2);
`output redirected...`> print(); # input placeholder
> junk := Multiply(capA-lambda*Iden, kv);
> %;
> righteigenvector := solve(subs(lambda = I*omega, {seq(junk(nn, 1), nn = 2 .. G+2)}), {seq(v[nn], nn = 2 .. G+2)});
`output redirected...`> print(); # input placeholder
> mo := `<,>`(1, seq(m[nn], nn = 2 .. G+2));
`output redirected...`> print(); # input placeholder
> junke := Multiply(capA-lambda*Iden, mo);
> righteigenvectorminus := solve(subs(lambda = -I*omega, {seq(junke(nn, 1), nn = 2 .. G+2)}), {seq(m[nn], nn = 2 .. G+2)});
`output redirected...`> print(); # input placeholder
> kw := `<|>`(1, Vector[row]([seq(w[nn], nn = 2 .. G+2)]));
`output redirected...`> print(); # input placeholder
> junk1 := Multiply(kw, capA-lambda*Iden);
> lefteigenvector := solve(subs(lambda = I*omega, {seq(junk1(1, nn), nn = 1 .. G+1)}), {seq(w[nn], nn = 2 .. G+2)});
`output redirected...`> print(); # input placeholder
> kwo := `<|>`(1, Vector[row]([seq(wo[nn], nn = 2 .. G+2)]));
`output redirected...`> print(); # input placeholder
> junke1 := Multiply(kwo, capA-lambda*Iden);
> lefteigenvectorminus := solve(subs(lambda = -I*omega, {seq(junke1(1, nn), nn = 1 .. G+1)}), {seq(wo[nn], nn = 2 .. G+2)});
`output redirected...`> print(); # input placeholder
> alpha := `<,>`(seq(diff(y[n](t), t), n = 1 .. G+2));
`output redirected...`> print(); # input placeholder
> beta := `<,>`(seq(y[n](t), n = 1 .. G+2));
`output redirected...`> print(); # input placeholder
> ee := alpha-Multiply(capB, beta);
> de[1] := ee(1, 1)-epsilon*h2*y[1](t)^2-epsilon^2*(h3*y[1](t)^3+y[2](t)*Delta);
`output redirected...`> print(); # input placeholder
> de[2] := ee(2, 1);
`output redirected...`> print(); # input placeholder
> for j to G do de[j+2] := ee(j+2, 1) end do;
`output redirected...`> print(); # input placeholder
> vm11 := Vector(G+2, [1, seq(v[nn], nn = 2 .. G+2)]);
`output redirected...`> print(); # input placeholder
> vm12 := Vector(G+2, [1, seq(m[nn], nn = 2 .. G+2)]);
`output redirected...`> print(); # input placeholder
> wm11 := Vector(G+2, [1, seq(w[nn], nn = 2 .. G+2)]);
# 
# 
> wm12 := Vector(G+2, [1, seq(wo[nn], nn = 2 .. G+2)]);
> for j to G+2 do for i to G+2 do expand(subs(y[i](t) = y[i](seq(T[k](t), k = 0 .. 2)), de[j])); expand(subs({seq(diff(T[k](t), t) = epsilon^k, k = 0 .. 2)}, %)); convert(taylor(expand(subs({seq(T[k](t) = T[k], k = 0 .. 2)}, %)), epsilon = 0, 3), polynom); de[j] := subs((D[1](y[j]))(seq(T[mmm], mmm = 0 .. 2)) = diff(y[j](seq(T[s], s = 0 .. 2)), T[0]), (D[2](y[j]))(seq(T[mmm], mmm = 0 .. 2)) = diff(y[j](seq(T[s], s = 0 .. 2)), T[1]), (D[3](y[j]))(seq(T[mmm], mmm = 0 .. 2)) = diff(y[j](seq(T[s], s = 0 .. 2)), T[2]), %) end do end do;
> for ff to G+2 do for j to G+2 do de[ff] := convert(taylor(expand(subs(y[j](seq(T[ss], ss = 0 .. 2)) = add(y[j, tt](T[0], T[1], T[2])*epsilon^tt, tt = 0 .. 2), de[ff])), epsilon = 0, 3), polynom) end do end do;
`output redirected...`> print(); # input placeholder
> for ff to G+2 do for j to G+2 do de[ff] := convert(taylor(combine(expand(subs(y[j, 0](seq(T[ss], ss = 0 .. 2)) = B1(T[1], T[2])*vm11(j, 1)*exp(I*omega*T[0])+B2(T[1], T[2])*vm12(j, 1)*exp(-I*omega*T[0]), de[ff]))), epsilon, 3), polynom) end do end do;
> for ff to G+2 do h1[ff] := coeff(coeff(de[ff], epsilon, 1), exp((1.*I)*omega*T[0])) end do;
> for ff to G+2 do h[ff] := coeff(coeff(de[ff], epsilon, 1), exp(-(1.*I)*omega*T[0])) end do;
> junk1 := sum('wm11[kk]*h1[kk]', 'kk' = 1 .. G+2);
> junk2 := sum('wm12[kk]*h[kk]', 'kk' = 1 .. G+2);
> eqs := solve({junk1, junk2}, {diff(B1(T[1], T[2]), T[1]), diff(B2(T[1], T[2]), T[1])});
??> print(); # input placeholder
> for ff to G+2 do t1[ff] := coeff(de[ff], epsilon, 1) end do;
> for j to G+2 do z1sol[j] := {y[j, 1](T[0], T[1], T[2]) = B1(T[1], T[2])^2*s[3*j-2]*exp((2*I)*omega*T[0])+B2(T[1], T[2])^2*s[3*j-1]*exp(-(2*I)*omega*T[0])+s[3*j]*B1(T[1], T[2])*B2(T[1], T[2])} end do;
> for ff to G+2 do t1[ff] := subs(eqs, convert(taylor(combine(expand(subs(seq(z1sol[j], j = 1 .. G+2), t1[ff]))), epsilon, 3), polynom)) end do;
> for ff to G+2 do t1[ff] := subs(exp((2*I)*omega*T[0]) = ev2, exp((2.*I)*omega*T[0]) = ev2, exp(-(2*I)*omega*T[0]) = en2, exp(-(2.*I)*omega*T[0]) = en2, exp((1.*I)*omega*T[0]) = ev1, exp(I*omega*T[0]) = ev1, exp(-(1.*I)*omega*T[0]) = en1, exp(-I*omega*T[0]) = en1, t1[ff]) end do;
> for rr to G+2 do eqe(3*rr-2) := coeff(t1[rr], ev2); eqe(3*rr-1) := coeff(t1[rr], en2); eqe(3*rr) := subs(ev2 = 0, en2 = 0, ev1 = 0, en1 = 0, t1[rr]) end do;
> sol := solve({seq(eqe(mm), mm = 1 .. 3*(G+2))}, {seq(s[nn], nn = 1 .. 3*(G+2))});
??> print(); # input placeholder
> for ff to G+2 do t2[ff] := simplify(subs(seq(z1sol[j], j = 1 .. G+2), coeff(de[ff], epsilon, 2))) end do; for ff to G+2 do t2[ff] := subs(exp((2*I)*omega*T[0]) = ev2, exp((2.*I)*omega*T[0]) = ev2, exp(-(2*I)*omega*T[0]) = en2, exp(-(2.*I)*omega*T[0]) = en2, exp((1.*I)*omega*T[0]) = ev1, exp(I*omega*T[0]) = ev1, exp(-(1.*I)*omega*T[0]) = en1, exp(-I*omega*T[0]) = en1, t2[ff]) end do;
> for ff to G+2 do h1[ff] := coeff(t2[ff], ev1) end do;
> for ff to G+2 do h[ff] := coeff(t2[ff], en1) end do;
> W1 := subs(sol, ee2, righteigenvector, righteigenvectorminus, lefteigenvector, lefteigenvectorminus, sum('wm11[kk]*h1[kk]', 'kk' = 1 .. G+2));
> W2 := subs(sol, ee2, righteigenvector, righteigenvectorminus, lefteigenvector, lefteigenvectorminus, sum('wm12[kk]*h[kk]', 'kk' = 1 .. G+2));
> junk11 := subs(phi(T[1], T[2]) = phi, R(T[1], T[2]) = R, subs(diff(R(T[1], T[2]), T[2]) = Rdot, diff(phi(T[1], T[2]), T[2]) = phidot, eval(subs(B1(T[1], T[2]) = (1/2)*R(T[1], T[2])*exp((1.*I)*phi(T[1], T[2])), B2(T[1], T[2]) = (1/2)*R(T[1], T[2])*exp(-(1.*I)*phi(T[1], T[2])), W1))));
> junk22 := subs(phi(T[1], T[2]) = phi, R(T[1], T[2]) = R, subs(diff(R(T[1], T[2]), T[2]) = Rdot, diff(phi(T[1], T[2]), T[2]) = phidot, eval(subs(B1(T[1], T[2]) = (1/2)*R(T[1], T[2])*exp((1.*I)*phi(T[1], T[2])), B2(T[1], T[2]) = (1/2)*R(T[1], T[2])*exp(-(1.*I)*phi(T[1], T[2])), W2))));
> eqsubs1 := {CR11 = coeff(subs(Rdot = 0, phidot = 0, junk11), R), CR13 = coeff(subs(Rdot = 0, phidot = 0, junk11), R^3), Cpd1 = simplify(coeff(junk11, phidot)), Crd1 = simplify(coeff(junk11, Rdot))};
> eqsubs2 := {CR21 = coeff(subs(Rdot = 0, phidot = 0, junk22), R), CR23 = coeff(subs(Rdot = 0, phidot = 0, junk22), R^3), Cpd2 = simplify(coeff(junk22, phidot)), Crd2 = simplify(coeff(junk22, Rdot))};
> eqsf := {Crd1*Rdot+Cpd1*phidot+CR11*R+CR13*R^3, Crd2*Rdot+Cpd2*phidot+CR21*R+CR23*R^3};
> junksol := solve(eqsf, {Rdot, phidot});
> Rsol := subs(junksol, Rdot);
> phisol := subs(junksol, phidot);
> Rsol := subs(eqsubs1, eqsubs2, Rsol);
> %;
> phisol := subs(eqsubs1, eqsubs2, phisol);
> tyu := expand(simplify(subs(omega = 1, h1 = 0.646e-1, h2 = 0.293e-1, h3 = -0.98e-1, Rsol)));
??> print(); # input placeholder
> limcye := evalc(Re(tyu));
??> print(); # input placeholder
> limcy := solve(limcye, R);
??> print(); # input placeholder
>