Try this —
syms yApprox(X1) a0 a1 a2 a3 a4 a5 a6
syms a2 a3 a4 a5 a6
EQ_01=a2 +9*a3 + 72*a4 + 540*a5 + 3888*a6 == -121/36000;
EQ_02=a3+ a2/12 +9*a4 +72*a5 +540*a6 == -1049/7776000;
EQ_03=a4 + (5*a2)/648 + (5*a3)/48 + (25*a5)/3 + (450*a6)/7 == -2987/419904000;
EQ_04 =a5 + (7*a2)/8640 + (7*a3)/600 + (7*a4)/60 + (63*a6)/8 == -19789/41990400000;
EQ_05 =a6 + a2/10800 + a3/720 + a4/70 + a5/8 == -24163/661348800000;
sol = solve([EQ_01, EQ_02, EQ_03, EQ_04, EQ_05],[a2,a3,a4,a5,a6])
a2 = vpa(sol.a2)
a3 = vpa(sol.a3)
a4 = vpa(sol.a4)
a5 = vpa(sol.a5)
a6 = vpa(sol.a6)
yApprox = a0 + a1*X1 + a2*X1^2 + a3*X1^3+ a4*X1^4 + a5*X1^5 + a6*X1^6;
yApprox = subs(yApprox)
yApprox = vpa(yApprox, 6) % Easier To Read, Retains Full Precison
.
