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.
