syms x a b c
f = x+(1/2).*a.*x.^2+0.3E0.*x.^3+(1/24).*(2.*a+(-0.1E0).*b+(-0.1E0).*c).*x.^4+(1/60).*(0.9E0+(1/2).*a.^2).*x.^5+(1/120).*( ...
0.6E0.*a+0.118333E-1.*b+(1/2).*(2.*a+(-0.1E0).*b+(-0.1E0).*c)+(1/2).*((-2).*a+0.1E0.*b+0.1E0.*c)+0.366667E-2.*c).*x.^6;
g = 1+b.*x+(-0.118333E0).*b.*x.^3+(-0.295833E-1).*a.*b.*x.^4+0.19525E-2.*b.*x.^5+(-0.236667E-1).*((-0.295833E0).*a.*b+ ...
(1/24).*b.*(2.*a+(-0.1E0).*b+(-0.1E0).*c)).*x.^6;
h = 1+c.*x+(-0.366667E-1).*c.*x.^3+(-0.916667E-2).*a.*c.*x.^4+(-0.209E-2).*c.*x.^5+(-0.733333E-2).*(( ...
-0.916667E-1).*a.*c+(1/24).*(2.*a+(-0.1E0).*b+(-0.1E0).*c).*c).* x.^6;
f22 = pade(diff(f),x,0,'Order',[2 2]);
g22 = pade(g,x,0,'Order',[2 2]);
h22 = pade(h,x,0,'Order',[2 2]);
f1 = limit(f22,x,Inf);
g1 = limit(g22,x,Inf);
h1 = limit(h22,x,Inf);
eqns = [f1 == 0, g1 == 0, h1 == 0];
S2 = vpasolve(eqns,[a b c]);
a_val = S2.a;
b_val = S2.b;
c_val = S2.c;
vars = [a b c];
a_val = a_val(real(a_val)&imag(a_val) == 0);
b_val = b_val(real(b_val)&imag(b_val) == 0);
c_val = c_val(real(c_val)&imag(c_val) == 0);
F = vpa (subs ([f, g, h], [a, b, c], [a_val(1) b_val(1) c_val(1)]) );
f = F(1); g = F(2); h = F(3);
disp([f g h])
Best Answer