About solving a system of equations by 'For' loop

I've written 4 equations within a for loop for 10 steps. But there's the following error:

Attempted to access C(1); index out of bounds because numel(C)=0.

Error in ROW (line 59)

k3=[C(1);C(2);C(3);C(4)];

When I use "double()" for any equation, the answer is not satisfactory for me. I was wondering if you could help me about my problem.

In the following, J1,J2,J3,J4, f1,f2,f3,f4 are knowns.

My currenct script looks like this:

h=1/1000;
y0=[2;2;1;0];
z0=10;
for i=1:10
   syms k11 k12 k13 k14 l11
k1=[k11;k12;k13;k14];
eqns = [h*(f1(y0(1),y0(2),y0(3),y0(4),z0)+0.44*(J1(1,:)*k1 +J2(1,:)*l11))-k11==0,
   h*(f2(y0(1),y0(2),y0(3),y0(4),z0)+0.44*(J1(2,:)*k1 +J2(2,:)*l11))-k12==0,
    h*(f3(y0(1),y0(2),y0(3),y0(4),z0)+0.44*(J1(3,:)*k1 +J2(3,:)*l11))-k13 ==0 ,
  h*(f4(y0(1),y0(2),y0(3),y0(4),z0)+0.44*(J1(4,:)*k1 +J2(4,:)*l11))-k14 ==0,
   g(y0(1),y0(2),y0(3),y0(4),z0)+0.44*(J3*k1 +J4*l11)==0];
vars = [k11, k12, k13, k14,  l11];
[solk11,solk12,solk13,solk14,soll11] = solve(eqns,vars);
A =double( [solk11,solk12,solk13,solk14,soll11]);
k1=[A(1);A(2);A(3);A(4)];
l11=A(5);


 syms k21 k22 k23 k24 l21
k2=[k21;k22;k23;k24];
eqns = [h*(f1(y0(1)+E(5)*k1(1),y0(2)+E(5)*k1(2),y0(3)+E(5)*k1(3),y0(4)+E(5)*k1(4),z0+E(5)*l11)+E(11)*(J1(1,:)*k1 +J2(1,:)*l11)+0.44*(J1(1,:)*k2 +J2(1,:)*l21))-k21==0,
h*(f2(y0(1)+E(5)*k1(1),y0(2)+E(5)*k1(2),y0(3)+E(5)*k1(3),y0(4)+E(5)*k1(4),z0+E(5)*l11)+E(11)*(J1(2,:)*k1 +J2(2,:)*l11)+0.44*(J1(2,:)*k2 +J2(2,:)*l21))-k22==0,
h*(f3(y0(1)+E(5)*k1(1),y0(2)+E(5)*k1(2),y0(3)+E(5)*k1(3),y0(4)+E(5)*k1(4),z0+E(5)*l11)+E(11)*(J1(3,:)*k1 +J2(3,:)*l11)+0.44*(J1(3,:)*k2 +J2(3,:)*l21))-k23 ==0 ,
h*(f4(y0(1)+E(5)*k1(1),y0(2)+E(5)*k1(2),y0(3)+E(5)*k1(3),y0(4)+E(5)*k1(4),z0+E(5)*l11)+E(11)*(J1(4,:)*k1 +J2(4,:)*l11)+0.44*(J1(4,:)*k2 +J2(4,:)*l21))-k24 ==0,
g(y0(1)+E(5)*k1(1),y0(2)+E(5)*k1(2),y0(3)+E(5)*k1(3),y0(4)+E(5)*k1(4),z0+E(5)*l11)+E(11)*(J3*k1 +J4*l11)+0.44*(J3*k2 +J4*l21)==0];
vars = [k21, k22, k23, k24,  l21];
[solk21,solk22,solk23,solk24,soll21] = solve(eqns,vars);
B =double( [solk21,solk22,solk23,solk24,soll21]);
k2=[B(1);B(2);B(3);B(4)];
l21=B(5);



 syms k31 k32 k33 k34 l31
k3=[k31;k32;k33;k34];
eqns = [h*(f1(y0(1)+E(6)*k1(1)+E(7)*k2(1),y0(2)+E(6)*k1(2)+E(7)*k2(2),y0(3)+E(6)*k1(3)+E(7)*k2(3),y0(4)+E(6)*k1(4)+E(7)*k2(4),z0+E(6)*l11+E(7)*l21)+E(12)*(J1(1,:)*k1 +J2(1,:)*l11)+E(13)*(J1(1,:)*k2 +J2(1,:)*l21)+0.44*(J1(1,:)*k3 +J2(1,:)*l31))-k31==0,
h*(f2(y0(1)+E(6)*k1(1)+E(7)*k2(1),y0(2)+E(6)*k1(2)+E(7)*k2(2),y0(3)+E(6)*k1(3)+E(7)*k2(3),y0(4)+E(6)*k1(4)+E(7)*k2(4),z0+E(6)*l11+E(7)*l21)+E(12)*(J1(2,:)*k1 +J2(2,:)*l11)+E(13)*(J1(2,:)*k2 +J2(2,:)*l21)+0.44*(J1(2,:)*k3 +J2(2,:)*l31))-k32==0,
h*(f3(y0(1)+E(6)*k1(1)+E(7)*k2(1),y0(2)+E(6)*k1(2)+E(7)*k2(2),y0(3)+E(6)*k1(3)+E(7)*k2(3),y0(4)+E(6)*k1(4)+E(7)*k2(4),z0+E(6)*l11+E(7)*l21)+E(12)*(J1(3,:)*k1 +J2(3,:)*l11)+E(13)*(J1(3,:)*k2 +J2(3,:)*l21)+0.44*(J1(3,:)*k3 +J2(3,:)*l31))-k33==0 ,
h*(f4(y0(1)+E(6)*k1(1)+E(7)*k2(1),y0(2)+E(6)*k1(2)+E(7)*k2(2),y0(3)+E(6)*k1(3)+E(7)*k2(3),y0(4)+E(6)*k1(4)+E(7)*k2(4),z0+E(6)*l11+E(7)*l21)+E(12)*(J1(4,:)*k1 +J2(4,:)*l11)+E(13)*(J1(4,:)*k2 +J2(4,:)*l21)+0.44*(J1(4,:)*k3 +J2(4,:)*l31))-k34 ==0,
g(y0(1)+E(6)*k1(1)+E(7)*k2(1),y0(2)+E(6)*k1(2)+E(7)*k2(2),y0(3)+E(6)*k1(3)+E(7)*k2(3),y0(4)+E(6)*k1(4)+E(7)*k2(4),z0+E(6)*l11+E(7)*l21)+E(11)*(J3*k1 +J4*l11)+E(13)*(J3*k2 +J4*l21)+0.44*(J3*k3 +J4*l31)==0];
vars = [k31, k32, k33, k34,  l31];
[solk31,solk32,solk33,solk34,soll31] = solve(eqns,vars);
C =double([solk31,solk32,solk33,solk34,soll31]);
k3=[C(1);C(2);C(3);C(4)];
l31=C(5);


syms k41 k42 k43 k44 l41
k4=[k41;k42;k43;k44];
eqns = [h*(f1(y0(1)+E(8)*k1(1)+E(8)*k2(1)+E(10)*k3(1),y0(2)+E(8)*k1(2)+E(8)*k2(2)+E(10)*k3(2),y0(3)+E(8)*k1(3)+E(8)*k2(3)+E(10)*k3(3),y0(4)+E(8)*k1(4)+E(8)*k2(4)+E(10)*k3(4),z0+E(8)*l11+E(9)*l21+E(10)*l31)+E(14)*(J1(1,:)*k1 +J2(1,:)*l11)+E(15)*(J1(1,:)*k2 +J2(1,:)*l21)+E(16)*(J1(1,:)*k3 +J2(1,:)*l31)+0.44*(J1(1,:)*k4 +J2(1,:)*l41))-k41==0,
h*(f2(y0(1)+E(8)*k1(1)+E(8)*k2(1)+E(10)*k3(1),y0(2)+E(8)*k1(2)+E(8)*k2(2)+E(10)*k3(2),y0(3)+E(8)*k1(3)+E(8)*k2(3)+E(10)*k3(3),y0(4)+E(8)*k1(4)+E(8)*k2(4)+E(10)*k3(4),z0+E(8)*l11+E(9)*l21+E(10)*l31)+E(14)*(J1(2,:)*k1 +J2(2,:)*l11)+E(15)*(J1(2,:)*k2 +J2(2,:)*l21)+E(16)*(J1(2,:)*k2 +J2(2,:)*l21)+0.44*(J1(2,:)*k4 +J2(2,:)*l41))-k42==0,
h*(f3(y0(1)+E(8)*k1(1)+E(8)*k2(1)+E(10)*k3(1),y0(2)+E(8)*k1(2)+E(8)*k2(2)+E(10)*k3(2),y0(3)+E(8)*k1(3)+E(8)*k2(3)+E(10)*k3(3),y0(4)+E(8)*k1(4)+E(8)*k2(4)+E(10)*k3(4),z0+E(8)*l11+E(9)*l21+E(10)*l31)+E(14)*(J1(3,:)*k1 +J2(3,:)*l11)+E(15)*(J1(3,:)*k2 +J2(3,:)*l21)+E(16)*(J1(3,:)*k3 +J2(3,:)*l31)+0.44*(J1(3,:)*k4 +J2(3,:)*l41))-k43 ==0 ,
h*(f4(y0(1)+E(8)*k1(1)+E(8)*k2(1)+E(10)*k3(1),y0(2)+E(8)*k1(2)+E(8)*k2(2)+E(10)*k3(2),y0(3)+E(8)*k1(3)+E(8)*k2(3)+E(10)*k3(3),y0(4)+E(8)*k1(4)+E(8)*k2(4)+E(10)*k3(4),z0+E(8)*l11+E(9)*l21+E(10)*l31)+E(14)*(J1(4,:)*k1 +J2(4,:)*l11)+E(15)*(J1(4,:)*k2 +J2(4,:)*l21)+E(16)*(J1(4,:)*k3 +J2(4,:)*l31)+0.44*(J1(4,:)*k4 +J2(4,:)*l41))-k44 ==0,
g(y0(1)+E(8)*k1(1)+E(8)*k2(1)+E(10)*k3(1),y0(2)+E(8)*k1(2)+E(8)*k2(2)+E(10)*k3(2),y0(3)+E(8)*k1(3)+E(8)*k2(3)+E(10)*k3(3),y0(4)+E(8)*k1(4)+E(8)*k2(4)+E(10)*k3(4),z0+E(8)*l11+E(9)*l21+E(10)*l31)+E(14)*(J3*k1 +J4*l11)+E(15)*(J3*k2 +J4*l21)+E(16)*(J3*k3 +J4*l31)+0.44*(J3*k4 +J4*l41)==0];
vars = [k41, k42, k43, k44,  l41];
[solk41,solk42,solk43,solk44,soll41] = solve(eqns,vars);
D =double( [solk41,solk42,solk43,solk44,soll41]);
k4=[D(1);D(2);D(3);D(4)];
l41=D(5);
y0=y0+E(1)*k1+E(2)*k2+E(3)*k3+E(4)*k4;
z0=z0+E(1)*l11+E(2)*l21+E(3)*l31+E(4)*l41;
end

[Walter Roberson]

How many columns do the J variables have?

[Mojtaba Mohareri]

J1, J2,J3,J4 are 4*4, 4*1, 1*4* ,1*1 matrices, respectively.

[Walter Roberson]

It appears that I will need f1 function, and J* values, and E values, and possibly some others, in order to run the code to test.

[Mojtaba Mohareri]

Thank you for your help. Yes, of course.

E=[0.3285609536316354, -0.5785609536316354, 0.25, 1, 2.218787467653286, 0, 0, 1.208587690772214, 0.07511610241919324, 0.5, -2.218787467653286, -0.09461966143940745, -0.007913526735718213, -1.870323744195384, -0.09624340112825115 , 0.2726301276675511  ];
f1(y1,y2,y3,y4,z1)=(-(z1^3)/(y3^2))*(3*(y2-y1+y3^(-1)-z1/10)^2+(1/5)*(y2-y1+y3^(-1)-(z1)/10))-y4;
f2(y1,y2,y3,y4,z1)=1/10*z1-y4;
f3(y1,y2,y3,y4,z1)=(z1^3)*(3*(y2-y1+y3^(-1)-z1/10)^2+(1/5)*(y2-y1+y3^(-1)-(z1)/10));
f4(y1,y2,y3,y4,z1)=y1-y3^(-1);
g(y1,y2,y3,y4,z1)=(y1-y3^(-1))^2+y4^2-1/10*z1;

J1(y1,y2,y3,y4,z1)=jacobian([f1,f2,f3,f4],[y1,y2,y3,y4]);
J1=J1(2,2,1,0,10);
J2(y1,y2,y3,y4,z1)=jacobian([f1,f2,f3,f4],[z1]);
J2=J2(2,2,1,0,10);
J3(y1,y2,y3,y4,z1)=jacobian([g],[y1,y2,y3,y4]);
J3=J3(2,2,1,0,10);
J4(y1,y2,y3,y4,z1)=jacobian([g],[z1]);
J4=J4(2,2,1,0,10);
y0=[2;2;1;0];
z0=10;

Answer 1

[Walter Roberson]

J1(1,:) is a vector of 4 elements. J2(1,:) is a scalar. When you use both in a single == expression then you are defining 4 simultaneous equations that vary only in J1 values, but you expect that single values of each of the variables will be able to satisfy all of those 4 equations simultaneously. That can only work if the defined equations are degenerate, where the J1 term is being multiplied by 0.

[Mojtaba Mohareri]

Could you explain to me what to do? I don't really understand it.

[Walter Roberson]

The first expression of your equations matrix uses J1(1,:) and you say that J1 is a 4x4 array. J1(1,:) would expand to a 1x4 vector. Your first expression is vector valued and so defines 4 equations all of which need to hold at the same time. Is that your intention?

[Mojtaba Mohareri]

J1 is a row vector, for example J1=[1,1,1,1] and k1 is also a column vector [k11;k12;k13;k14]. So J1(1,:)*k1 is k11+k12+k13+k14 and is not vector valued.

[Walter Roberson]

No, J1 is 4*4. You told us that when I asked you about the sizes. https://www.mathworks.com/matlabcentral/answers/455017-about-solving-a-system-of-equations-by-for-loop#comment_691120

[Mojtaba Mohareri]

I'm so sorry, your are right. J1 a is 4*4 matrix, but J1(1,:) is a vector.

[Mojtaba Mohareri]

So, as I told you, J1(1,:)*k1 is a scaler.

[Mojtaba Mohareri]

Dear Walter Roberson. Can I solve my problem with other code? I want just to obtain y0 and z0 after 1000 steps and it's not important for me which code I use. Unfortunately I don't know matlab properly. Could you help me to get the answer?

[Walter Roberson]

It appears that the double() you are using for efficiency is leading to enough loss of precision that you are getting Inf solutions early on, which lead to NaN in further calculation, which lead to your eqns in the C block all immediately evaluating to FALSE before even reaching the solve(), which leads to the solve() returning empty, which leads to subscript out of range.

When I remove those double(), the calculations are much much much much slower, and iterating i to 1000 would probably not be practical.

I will see what I can come up with.

[Walter Roberson]

When I change the double() to vpa() with 16 digits (equivalent precision but larger range) then by i = 10, the A matrix is composed of values that have magnitude greater than 10^one million and the symbolic numeric solver starts giving up. That is, your equations are overflowing, and it is not due to round-off (same problem happens with 32 digits for example.)

Everything looks fine in the A matrix until i = 6, at which point you start seeing values > 10 after having been no more than +/- 2 before that. i = 7 gives you A matrix values about 1E65. i = 8 gives you A matrix values about 1E8000 . If you were working in double precision, that would be the point at which the values would overflow to infinity and everything after that would start to fail.

You will need to go back and re-check the equations. With the way they are now, you cannot get through 10 iterations without overflowing symbolic numeric range, so there is no way you could get through 1000 iterations.

Your f1 involves z1^4 and the values that result from f1 generally feed back to become the arguments to f1. When z1 becomes greater than 1, that leads to a chain of 4th powers in increases in range, leading to rapid overflow.

[Mojtaba Mohareri]

Thak you for your exact explanation and your help. As you told me, I get the answer by i=6 but not for values i greater than 6. Yes, I'll re-check my equations. I think I can remove somehow z1 from the equations but I have to check them again.

Thank you very much for taking the trouble to help me. I really appreciate your help.

[Mojtaba Mohareri]

Dear Walter Roberson, I simplify my problem by assuming z0 to be constant 10 during For Loop. I've also changed the double() to vpa() as you told. Now, I get the answer by i=200, but for i=1000. there's an error as

Error using subsref

Index exceeds matrix dimensions. Error in sym/subsref (line 771)

R_tilde = builtin('subsref',L_tilde,Idx);

Error in ROW (line 46)

k2=[B(1);B(2);B(3);B(4)];

Could you tell me what to do for it?

[Walter Roberson]

Please post your current code. When I try your code as posted before except with the change to z0 commented out, then it stops at i = 102 with "integer too large in context" because it has reached values that exceed 10 to the three and a half million power.

[Mojtaba Mohareri]

Yes, you are right, but I don't know what the problem is. In your opinion, Is there something wrong with my equations? or with my codes? I re-check my equations every time:(

Yes, sure. My current code:

syms b1 b2 b_3 b_4  a21 a31 a32 a41 a42 a43 g21 g31 g32 g41 g42 g43  w12 w22 y1 y2 y3 y4 z1 k11 k12 k13 k14 k21 k22 k23 k24 k31 k32 k33 k34 k41 k42 k43 k44 l11 l21 l31 l41 x
k=[k11;k12;k13;k14];
a=0;
h=1/1000;
E=[0.3285609536316354, -0.5785609536316354, 0.25, 1, 2.218787467653286, 0, 0, 1.208587690772214, 0.07511610241919324, 0.5, -2.218787467653286, -0.09461966143940745, -0.007913526735718213, -1.870323744195384, -0.09624340112825115 , 0.2726301276675511  ];
f1(y1,y2,y3,y4,z1)=(-(z1^3)/(y3^2))*(3*(y2-y1+y3^(-1)-z1/10)^2+(1/5)*(y2-y1+y3^(-1)-(z1)/10))-y4;
f2(y1,y2,y3,y4,z1)=1/10*z1-y4;
f3(y1,y2,y3,y4,z1)=(z1^3)*(3*(y2-y1+y3^(-1)-(z1)/10)^2+(1/5)*(y2-y1+y3^(-1)-(z1)/10));
f4(y1,y2,y3,y4,z1)=y1-y3^(-1);
g(y1,y2,y3,y4,z1)=(y1-y3^(-1))^2+y4^2-1/10*z1;

J1(y1,y2,y3,y4,z1)=jacobian([f1,f2,f3,f4],[y1,y2,y3,y4]);
J1=J1(2,2,1,0,10);
J2(y1,y2,y3,y4,z1)=jacobian([f1,f2,f3,f4],[z1]);
J2=J2(2,2,1,0,10);
J3(y1,y2,y3,y4,z1)=jacobian([g],[y1,y2,y3,y4]);
J3=J3(2,2,1,0,10);
J4=-1/10;
y0=[2;2;1;0];
z0=10;


for i=1:500
   syms k11 k12 k13 k14 l11
k1=[k11;k12;k13;k14];
eqns = [h*(f1(y0(1),y0(2),y0(3),y0(4),z0)+0.44*(J1(1,:)*k1 +J2(1,:)*l11))-k11==0,
   h*(f2(y0(1),y0(2),y0(3),y0(4),z0)+0.44*(J1(2,:)*k1 +J2(2,:)*l11))-k12==0,
     h*(f3(y0(1),y0(2),y0(3),y0(4),z0)+0.44*(J1(3,:)*k1 +J2(3,:)*l11))-k13 ==0 ,
     h*(f4(y0(1),y0(2),y0(3),y0(4),z0)+0.44*(J1(4,:)*k1 +J2(4,:)*l11))-k14 ==0,
     g(y0(1),y0(2),y0(3),y0(4),z0)+0.44*(J3*k1 +J4*l11)==0];
vars = [k11, k12, k13, k14,  l11];
[solk11,solk12,solk13,solk14,soll11] = solve(eqns,vars);
A =vpa([solk11,solk12,solk13,solk14,soll11]);
k1=[A(1);A(2);A(3);A(4)];
l11=A(5);
 syms k21 k22 k23 k24 l21
k2=[k21;k22;k23;k24];
eqns = [h*(f1(y0(1)+E(5)*k1(1),y0(2)+E(5)*k1(2),y0(3)+E(5)*k1(3),y0(4)+E(5)*k1(4),z0)+E(11)*(J1(1,:)*k1 +J2(1,:)*l11)+0.44*(J1(1,:)*k2 +J2(1,:)*l21))-k21==0,
 h*(f2(y0(1)+E(5)*k1(1),y0(2)+E(5)*k1(2),y0(3)+E(5)*k1(3),y0(4)+E(5)*k1(4),z0)+E(11)*(J1(2,:)*k1 +J2(2,:)*l11)+0.44*(J1(2,:)*k2 +J2(2,:)*l21))-k22==0,
 h*(f3(y0(1)+E(5)*k1(1),y0(2)+E(5)*k1(2),y0(3)+E(5)*k1(3),y0(4)+E(5)*k1(4),z0)+E(11)*(J1(3,:)*k1 +J2(3,:)*l11)+0.44*(J1(3,:)*k2 +J2(3,:)*l21))-k23 ==0 ,
 h*(f4(y0(1)+E(5)*k1(1),y0(2)+E(5)*k1(2),y0(3)+E(5)*k1(3),y0(4)+E(5)*k1(4),z0)+E(11)*(J1(4,:)*k1 +J2(4,:)*l11)+0.44*(J1(4,:)*k2 +J2(4,:)*l21))-k24 ==0,
 g(y0(1)+E(5)*k1(1),y0(2)+E(5)*k1(2),y0(3)+E(5)*k1(3),y0(4)+E(5)*k1(4),z0)+E(11)*(J3*k1 +J4*l11)+0.44*(J3*k2 +J4*l21)==0];
vars = [k21, k22, k23, k24,  l21];
[solk21,solk22,solk23,solk24,soll21] = solve(eqns,vars);
B =vpa( [solk21,solk22,solk23,solk24,soll21]);
k2=[B(1);B(2);B(3);B(4)];
l21=B(5);
 syms k31 k32 k33 k34 l31
k3=[k31;k32;k33;k34];
eqns = [h*(f1(y0(1)+E(6)*k1(1)+E(7)*k2(1),y0(2)+E(6)*k1(2)+E(7)*k2(2),y0(3)+E(6)*k1(3)+E(7)*k2(3),y0(4)+E(6)*k1(4)+E(7)*k2(4),z0)+E(12)*(J1(1,:)*k1 +J2(1,:)*l11)+E(13)*(J1(1,:)*k2 +J2(1,:)*l21)+0.44*(J1(1,:)*k3 +J2(1,:)*l31))-k31==0,
h*(f2(y0(1)+E(6)*k1(1)+E(7)*k2(1),y0(2)+E(6)*k1(2)+E(7)*k2(2),y0(3)+E(6)*k1(3)+E(7)*k2(3),y0(4)+E(6)*k1(4)+E(7)*k2(4),z0)+E(12)*(J1(2,:)*k1 +J2(2,:)*l11)+E(13)*(J1(2,:)*k2 +J2(2,:)*l21)+0.44*(J1(2,:)*k3 +J2(2,:)*l31))-k32==0,
h*(f3(y0(1)+E(6)*k1(1)+E(7)*k2(1),y0(2)+E(6)*k1(2)+E(7)*k2(2),y0(3)+E(6)*k1(3)+E(7)*k2(3),y0(4)+E(6)*k1(4)+E(7)*k2(4),z0)+E(12)*(J1(3,:)*k1 +J2(3,:)*l11)+E(13)*(J1(3,:)*k2 +J2(3,:)*l21)+0.44*(J1(3,:)*k3 +J2(3,:)*l31))-k33 ==0 ,
h*(f4(y0(1)+E(6)*k1(1)+E(7)*k2(1),y0(2)+E(6)*k1(2)+E(7)*k2(2),y0(3)+E(6)*k1(3)+E(7)*k2(3),y0(4)+E(6)*k1(4)+E(7)*k2(4),z0)+E(12)*(J1(4,:)*k1 +J2(4,:)*l11)+E(13)*(J1(4,:)*k2 +J2(4,:)*l21)+0.44*(J1(4,:)*k3 +J2(4,:)*l31))-k34 ==0,
g(y0(1)+E(6)*k1(1)+E(7)*k2(1),y0(2)+E(6)*k1(2)+E(7)*k2(2),y0(3)+E(6)*k1(3)+E(7)*k2(3),y0(4)+E(6)*k1(4)+E(7)*k2(4),z0)+E(11)*(J3*k1 +J4*l11)+E(13)*(J3*k2 +J4*l21)+0.44*(J3*k3 +J4*l31)==0];
vars = [k31, k32, k33, k34,  l31];
[solk31,solk32,solk33,solk34,soll31] = solve(eqns,vars);
C =vpa([solk31,solk32,solk33,solk34,soll31]);
k3=[C(1);C(2);C(3);C(4)];
l31=C(5);
syms k41 k42 k43 k44 l41
k4=[k41;k42;k43;k44];
eqns = [h*(f1(y0(1)+E(8)*k1(1)+E(8)*k2(1)+E(10)*k3(1),y0(2)+E(8)*k1(2)+E(8)*k2(2)+E(10)*k3(2),y0(3)+E(8)*k1(3)+E(8)*k2(3)+E(10)*k3(3),y0(4)+E(8)*k1(4)+E(8)*k2(4)+E(10)*k3(4),z0)+E(14)*(J1(1,:)*k1 +J2(1,:)*l11)+E(15)*(J1(1,:)*k2 +J2(1,:)*l21)+E(16)*(J1(1,:)*k3 +J2(1,:)*l31)+0.44*(J1(1,:)*k4 +J2(1,:)*l41))-k41==0,
h*(f2(y0(1)+E(8)*k1(1)+E(8)*k2(1)+E(10)*k3(1),y0(2)+E(8)*k1(2)+E(8)*k2(2)+E(10)*k3(2),y0(3)+E(8)*k1(3)+E(8)*k2(3)+E(10)*k3(3),y0(4)+E(8)*k1(4)+E(8)*k2(4)+E(10)*k3(4),z0)+E(14)*(J1(2,:)*k1 +J2(2,:)*l11)+E(15)*(J1(2,:)*k2 +J2(2,:)*l21)+E(16)*(J1(2,:)*k2 +J2(2,:)*l21)+0.44*(J1(2,:)*k4 +J2(2,:)*l41))-k42==0,
h*(f3(y0(1)+E(8)*k1(1)+E(8)*k2(1)+E(10)*k3(1),y0(2)+E(8)*k1(2)+E(8)*k2(2)+E(10)*k3(2),y0(3)+E(8)*k1(3)+E(8)*k2(3)+E(10)*k3(3),y0(4)+E(8)*k1(4)+E(8)*k2(4)+E(10)*k3(4),z0)+E(14)*(J1(3,:)*k1 +J2(3,:)*l11)+E(15)*(J1(3,:)*k2 +J2(3,:)*l21)+E(16)*(J1(3,:)*k3 +J2(3,:)*l31)+0.44*(J1(3,:)*k4 +J2(3,:)*l41))-k43 ==0 ,
h*(f4(y0(1)+E(8)*k1(1)+E(8)*k2(1)+E(10)*k3(1),y0(2)+E(8)*k1(2)+E(8)*k2(2)+E(10)*k3(2),y0(3)+E(8)*k1(3)+E(8)*k2(3)+E(10)*k3(3),y0(4)+E(8)*k1(4)+E(8)*k2(4)+E(10)*k3(4),z0)+E(14)*(J1(4,:)*k1 +J2(4,:)*l11)+E(15)*(J1(4,:)*k2 +J2(4,:)*l21)+E(16)*(J1(4,:)*k3 +J2(4,:)*l31)+0.44*(J1(4,:)*k4 +J2(4,:)*l41))-k44 ==0,
g(y0(1)+E(8)*k1(1)+E(8)*k2(1)+E(10)*k3(1),y0(2)+E(8)*k1(2)+E(8)*k2(2)+E(10)*k3(2),y0(3)+E(8)*k1(3)+E(8)*k2(3)+E(10)*k3(3),y0(4)+E(8)*k1(4)+E(8)*k2(4)+E(10)*k3(4),z0)+E(14)*(J3*k1 +J4*l11)+E(15)*(J3*k2 +J4*l21)+E(16)*(J3*k3 +J4*l31)+0.44*(J3*k4 +J4*l41)==0];
vars = [k41, k42, k43, k44,  l41];
[solk41,solk42,solk43,solk44,soll41] = solve(eqns,vars);
D =vpa( [solk41,solk42,solk43,solk44,soll41]);
k4=[D(1);D(2);D(3);D(4)];
l41=D(5);
y0=y0+E(1)*k1+E(2)*k2+E(3)*k3+E(4)*k4;
y{i}=y0;
w(i)=y0(4);
z0=10;
end

Thank you very much.