MATLAB: When I run the program inaccurate results appear What is the problem

nonlinear

syms t y1 y2 y3 y4  yn11 yn1 yr yr11 yr111   yr v y(v) x
format longe
h=0.1
r=1.7686999343761197980595002942562   %1694/911=1.859495060373216e+00
x=0
y0=0
y01=0.25;
y011=0;
fn=-y01+y0.*y01.*y011
fnr=-yr11+yr.*yr11.*yr111
fn1=-yn1+y1.*yn1.*yn11
fn2=-y3+y2.*y3.*y4
eq1 =-y2-(r - 2)/r*y0+ 2/(r*(r - 1))*yr+ (2*r - 4)/(r - 1)*y1 -(h.^3.*(r^4 - 8*r^3 + 22*r^2 - 23*r + 6))/(120*r)*fn+ (h.^3.*(- r^2 + 2*r + 3))/(60*r)*fnr -(h.^3.*(- r^3 + 4*r^2 + 9*r - 26))/60*fn1+ (h.^3.*(- r^3 + 2*r + 1))/120*fn2
 
eq2= -h.*y3 -(r - 3)/r*y0+ 3/(r*(r - 1))*yr+ (r - 4)/(r - 1)*y1 -(h.^3.*(3*r^4 - 24*r^3 + 66*r^2 - 67*r + 12))/(240*r)*fn -(h.^3.*(r^4 - 5*r^3 + 5*r^2 + 5*r - 4))/(40*r*(r^2 - 3*r + 2))*fnr -(h.^3.*(- 3*r^4 + 15*r^3 + 15*r^2 - 137*r + 116))/(120*(r - 1))*fn1+ (h.^3.*(- r^4 + 2*r^3 + 2*r^2 + 3*r - 12))/(80*(r - 2))*fn2
 
eq3=-h.^2.*y4+2/r*y0+ 2/(r*(r - 1))*yr -2/(r - 1)*y1+ (h.^3.*(- r^4 + 8*r^3 - 22*r^2 + 18*r + 5))/(120*r)*fn -(h.^3.*(r^4 - 5*r^3 + 5*r^2 + 5*r + 5))/(60*r*(r^2 - 3*r + 2))*fnr -(h.^3.*(- r^4 + 5*r^3 + 5*r^2 - 75*r + 77))/(60*(r - 1))*fn1+ (h.^3.*(- r^4 + 2*r^3 + 2*r^2 + 42*r - 81))/(120*(r - 2))*fn2
 
eq4=-h.*yn1 -(r - 1)/r*y0+ 1/(r*(r - 1))*yr+ (r - 2)/(r - 1)*y1 -(h.^3.*(r^4 - 8*r^3 + 22*r^2 - 21*r + 6))/(240*r)*fn+ (h.^3.*(- r^2 + 2*r + 3))/(120*r)*fnr -(h.^3.*(- r^3 + 4*r^2 + 9*r - 8))/120*fn1 -(h.^3.*(r^3 - 2*r + 1))/240*fn2
 
eq5=-h.^2.*yn11+ 2/r*y0+ 2/(r*(r - 1))*yr -2/(r - 1)*y1 -(h.^3.*(r^4 - 8*r^3 + 22*r^2 - 28*r + 10))/(120*r)*fn -(h.^3.*(r^4 - 5*r^3 + 5*r^2 + 5*r - 10))/(60*r*(r^2 - 3*r + 2))*fnr -(h.^3.*(- r^4 + 5*r^3 + 5*r^2 - 35*r + 22))/(60*(r - 1))*fn1+ (h.^3.*(- r^4 + 2*r^3 + 2*r^2 - 8*r + 4))/(120*(r - 2))*fn2
eq6=-h.*yr11+ (r - 1)/r*y0+ (2*r - 1)/(r*(r - 1))*yr -r/(r - 1)*y1+ (h.^3.*(2*r^4 - 13*r^3 + 28*r^2 - 22*r + 5))/240*fn+ (h.^3.*(4*r^3 - 15*r^2 + 10*r + 5))/(120*(r - 2))*fnr+ (h.^3.*r*(- 2*r^3 + 7*r^2 + 2*r - 3))/120*fn1+ (h.^3.*r*(2*r^4 - 5*r^3 + 2*r^2 + 2*r - 1))/(240*(r - 2))*fn2
 
eq7=-h.^2.*yr111+ 2/r*y0+ 2/(r*(r - 1))*yr -2/(r - 1)*y1+ (h.^3.*(4*r^4 - 22*r^3 + 38*r^2 - 22*r + 5))/(120*r)*fn -(h.^3.*(- 14*r^4 + 55*r^3 - 55*r^2 + 5*r + 5))/(60*r*(r^2 - 3*r + 2))*fnr -(h.^3.*(4*r^4 - 15*r^3 + 5*r^2 + 5*r - 3))/(60*(r - 1))*fn1+ (h.^3.*(4*r^4 - 8*r^3 + 2*r^2 + 2*r - 1))/(120*(r - 2))*fn2
 
eq8=-h.*y01 -(r + 1)/r*y0 -1/(r*(r - 1))*yr+ r/(r - 1)*y1+ (h.^3.*(r^3 - 8*r^2 + 22*r - 5))/240*fn+ (h.^3.*(r^3 - 5*r^2 + 5*r + 5))/(120*(r^2 - 3*r + 2))*fnr+ (h.^3.*r*(- r^3 + 5*r^2 + 5*r - 3))/(120*(r - 1))*fn1 -(h.^3.*r*(- r^3 + 2*r^2 + 2*r - 1))/(240*(r - 2))*fn2
 
eq9 =-h.^2.*y011+ 2/r*y0+ 2/(r*(r - 1))*yr -2/(r - 1)*y1 -(h.^3.*(r^4 - 8*r^3 + 22*r^2 + 22*r - 5))/(120*r)*fn -(h.^3.*(r^4 - 5*r^3 + 5*r^2 + 5*r + 5))/(60*r*(r^2 - 3*r + 2))*fnr -(h.^3.*(- r^4 + 5*r^3 + 5*r^2 + 5*r - 3))/(60*(r - 1))*fn1+ (h.^3.*(- r^4 + 2*r^3 + 2*r^2 + 2*r - 1))/(120*(r - 2))*fn2
 [y1,y2,y3,y4,yn1,yn11,yr,yr11,yr111]=vpasolve([eq1,eq2,eq3,eq4,eq5,eq6,eq7,eq8,eq9])

Best Answer

The best approach to this is likely not the vpasolve function, since it takes forever.
I have no idea of these are the ‘accurate’ results you want (I have no idea what that means), however you will get them much more quickly: 
syms t y1 y2 y3 y4  yn11 yn1 yr yr11 yr111   yr v y(v) x
format longe
h=0.1;
r=1.7686999343761197980595002942562;   %1694/911=1.859495060373216e+00
x=0;
y0=0;
y01=0.25;
y011=0;
fn=-y01+y0.*y01.*y011;
fnr=-yr11+yr.*yr11.*yr111;
fn1=-yn1+y1.*yn1.*yn11;
fn2=-y3+y2.*y3.*y4;
eq1 =-y2-(r - 2)/r*y0+ 2/(r*(r - 1))*yr+ (2*r - 4)/(r - 1)*y1 -(h.^3.*(r^4 - 8*r^3 + 22*r^2 - 23*r + 6))/(120*r)*fn+ (h.^3.*(- r^2 + 2*r + 3))/(60*r)*fnr -(h.^3.*(- r^3 + 4*r^2 + 9*r - 26))/60*fn1+ (h.^3.*(- r^3 + 2*r + 1))/120*fn2;
eq1 = simplifyFraction(eq1);
eq1d = vpa(eq1, 5)
eq2= -h.*y3 -(r - 3)/r*y0+ 3/(r*(r - 1))*yr+ (r - 4)/(r - 1)*y1 -(h.^3.*(3*r^4 - 24*r^3 + 66*r^2 - 67*r + 12))/(240*r)*fn -(h.^3.*(r^4 - 5*r^3 + 5*r^2 + 5*r - 4))/(40*r*(r^2 - 3*r + 2))*fnr -(h.^3.*(- 3*r^4 + 15*r^3 + 15*r^2 - 137*r + 116))/(120*(r - 1))*fn1+ (h.^3.*(- r^4 + 2*r^3 + 2*r^2 + 3*r - 12))/(80*(r - 2))*fn2;
eq2 = simplifyFraction(eq2);
eq2d = vpa(eq2, 5)
eq3=-h.^2.*y4+2/r*y0+ 2/(r*(r - 1))*yr -2/(r - 1)*y1+ (h.^3.*(- r^4 + 8*r^3 - 22*r^2 + 18*r + 5))/(120*r)*fn -(h.^3.*(r^4 - 5*r^3 + 5*r^2 + 5*r + 5))/(60*r*(r^2 - 3*r + 2))*fnr -(h.^3.*(- r^4 + 5*r^3 + 5*r^2 - 75*r + 77))/(60*(r - 1))*fn1+ (h.^3.*(- r^4 + 2*r^3 + 2*r^2 + 42*r - 81))/(120*(r - 2))*fn2;
eq3 = simplifyFraction(eq3);
eq3d = vpa(eq3, 5)
eq4=-h.*yn1 -(r - 1)/r*y0+ 1/(r*(r - 1))*yr+ (r - 2)/(r - 1)*y1 -(h.^3.*(r^4 - 8*r^3 + 22*r^2 - 21*r + 6))/(240*r)*fn+ (h.^3.*(- r^2 + 2*r + 3))/(120*r)*fnr -(h.^3.*(- r^3 + 4*r^2 + 9*r - 8))/120*fn1 -(h.^3.*(r^3 - 2*r + 1))/240*fn2;
eq4 = simplifyFraction(eq4);
eq4d = vpa(eq4, 5)
eq5=-h.^2.*yn11+ 2/r*y0+ 2/(r*(r - 1))*yr -2/(r - 1)*y1 -(h.^3.*(r^4 - 8*r^3 + 22*r^2 - 28*r + 10))/(120*r)*fn -(h.^3.*(r^4 - 5*r^3 + 5*r^2 + 5*r - 10))/(60*r*(r^2 - 3*r + 2))*fnr -(h.^3.*(- r^4 + 5*r^3 + 5*r^2 - 35*r + 22))/(60*(r - 1))*fn1+ (h.^3.*(- r^4 + 2*r^3 + 2*r^2 - 8*r + 4))/(120*(r - 2))*fn2;
eq5 = simplifyFraction(eq5);
eq5d = vpa(eq5, 5)
eq6=-h.*yr11+ (r - 1)/r*y0+ (2*r - 1)/(r*(r - 1))*yr -r/(r - 1)*y1+ (h.^3.*(2*r^4 - 13*r^3 + 28*r^2 - 22*r + 5))/240*fn+ (h.^3.*(4*r^3 - 15*r^2 + 10*r + 5))/(120*(r - 2))*fnr+ (h.^3.*r*(- 2*r^3 + 7*r^2 + 2*r - 3))/120*fn1+ (h.^3.*r*(2*r^4 - 5*r^3 + 2*r^2 + 2*r - 1))/(240*(r - 2))*fn2;
eq6 = simplifyFraction(eq6);
eq6d = vpa(eq6, 5)
eq7=-h.^2.*yr111+ 2/r*y0+ 2/(r*(r - 1))*yr -2/(r - 1)*y1+ (h.^3.*(4*r^4 - 22*r^3 + 38*r^2 - 22*r + 5))/(120*r)*fn -(h.^3.*(- 14*r^4 + 55*r^3 - 55*r^2 + 5*r + 5))/(60*r*(r^2 - 3*r + 2))*fnr -(h.^3.*(4*r^4 - 15*r^3 + 5*r^2 + 5*r - 3))/(60*(r - 1))*fn1+ (h.^3.*(4*r^4 - 8*r^3 + 2*r^2 + 2*r - 1))/(120*(r - 2))*fn2;
eq7 = simplifyFraction(eq7);
eq7d = vpa(eq7, 5)
eq8=-h.*y01 -(r + 1)/r*y0 -1/(r*(r - 1))*yr+ r/(r - 1)*y1+ (h.^3.*(r^3 - 8*r^2 + 22*r - 5))/240*fn+ (h.^3.*(r^3 - 5*r^2 + 5*r + 5))/(120*(r^2 - 3*r + 2))*fnr+ (h.^3.*r*(- r^3 + 5*r^2 + 5*r - 3))/(120*(r - 1))*fn1 -(h.^3.*r*(- r^3 + 2*r^2 + 2*r - 1))/(240*(r - 2))*fn2;
eq8 = simplifyFraction(eq8);
eq8d = vpa(eq8, 5)
eq9 =-h.^2.*y011+ 2/r*y0+ 2/(r*(r - 1))*yr -2/(r - 1)*y1 -(h.^3.*(r^4 - 8*r^3 + 22*r^2 + 22*r - 5))/(120*r)*fn -(h.^3.*(r^4 - 5*r^3 + 5*r^2 + 5*r + 5))/(60*r*(r^2 - 3*r + 2))*fnr -(h.^3.*(- r^4 + 5*r^3 + 5*r^2 + 5*r - 3))/(60*(r - 1))*fn1+ (h.^3.*(- r^4 + 2*r^3 + 2*r^2 + 2*r - 1))/(120*(r - 2))*fn2;
eq9 = simplifyFraction(eq9);
eq9d = vpa(eq9, 5)
eqsfcn = matlabFunction([eq1,eq2,eq3,eq4,eq5,eq6,eq7,eq8,eq9], 'Vars',{[y1,y2,y3,y4,yn1,yn11,yr,yr11,yr111]})
eqsol = fsolve(eqsfcn, rand(1,9))
eqsolc = num2cell(eqsol);
[y1,y2,y3,y4,yn1,yn11,yr,yr11,yr111] = eqsolc{:}
I used the simplifyFraction function to do just that with your equations, then the fsolve function to do the actual calculations. These changes make it significantly more efficient.
The results I got are: 
y1 =
     2.495835158216658e-02
y2 =
     4.966725033036157e-02
y3 =
     2.450145921764798e-01
y4 =
    -4.970810349993776e-02
yn1 =
     2.487509121241510e-01
yn11 =
    -2.496352514349295e-02
yr =
     4.398727143758653e-02
yr11 =
     2.460985502335340e-01
yr111 =
    -4.401564490668385e-02
You must be certain your equations are written as you want them to be in order to get ‘accurate’ values. We have no control over that.
Experiment to get the results you want.
Related Question