how can improve my code performanece, it uses a syms vector, i tried to delete y=t= sym(zeros(1, m + 3)); and it's faster but the solution of the function its wrong

function [fi, fi1, Ji, Ji1] = funcion(L,r1,s1,t1,C0,C1,C2,C3,FS,b,D0,D1,D2,AstG,S,R,DgS,ra,n)
%definimos las condiciones iniciales
m = 40;
y = sym(zeros(1, m + 3));
t = sym(zeros(1, m + 3));
    y(1)=C0;
    y(2)=C1;
    t(1)=C2;
    t(2)=C3;
     
for k = 3:m
        y(k)=((1+FS)*(k+1)*t(k+2)+(2*R*k+2*R)*t(k+1)+(R^2*(k-1)+2*R^2)*t(k)-(2*R*k*(k+1)+2*R*(k+1))*y(k+2)+(-R^2*(k-1)*k+D0*S-b*S*n-2*R^2*k)*y(k+1)+(D1*S-2*b*S*n*R)*y(k)+(D2*S-b*S*n*R^2)*y(k-1))/((1+DgS)*(k+1)*(k+2));
        t(k)=(-(4*S*R*k*(k+1)+4*S*R*(k+1))*t(k+2)+(1+FS-b*S*ra-6*S*R^2*(k-1)*k-12*S*R^2*k)*t(k+1)+(-4*S*R^3*(k-2)*(k-1)+2*R-4*b*S*ra*R-12*S*R^3*(k-1))*t(k)+(-S*R^4*(k-3)*(k-2)+R^2-6*b*S*ra*R^2-4*S*R^4*(k-2))*t(k-1)-(4*b*S*ra*R^3)*t(k-2)-(b*S*ra*R^4)*t(k-3)-(1+FS)*(k+1)*y(k+2)-2*R*k*y(k+1)-R^2*(k-1)*y(k))/(S*(k+1)*(k+2));
end 
j = 0:m+2;
fi = sum(y);
fi1 = y * j';
Ji = sum(t);
Ji1 = t * j';
end       

Answer 1

[Torsten]

You can't use t(k), t(k+1), t(k+2), y(k), y(k+1) and y(k+2) to define y(k). The values for these y's and t's are not yet known.

Same for t(k).

Answer 2

[Walter Roberson]

        t(k)=(-(4*S*R*k*(k+1)+4*S*R*(k+1))*t(k+2)+(1+FS-b*S*ra-6*S*R^2*(k-1)*k-12*S*R^2*k)*t(k+1)+(-4*S*R^3*(k-2)*(k-1)+2*R-4*b*S*ra*R-12*S*R^3*(k-1))*t(k)+(-S*R^4*(k-3)*(k-2)+R^2-6*b*S*ra*R^2-4*S*R^4*(k-2))*t(k-1)-(4*b*S*ra*R^3)*t(k-2)-(b*S*ra*R^4)*t(k-3)-(1+FS)*(k+1)*y(k+2)-2*R*k*y(k+1)-R^2*(k-1)*y(k))/(S*(k+1)*(k+2));
       %                                                                                                                                                                                                                                                  ^^^^^^
       

I pointed to a place in the code where you use t(k-3) when k starts from 3.

[Andres Cordoba Ordoñez]

Actually is defined in my scrip, i show this like this to be brief

[Walter Roberson]

Your posted code:

for k = 3:m
        y(k)=((1+FS)*(k+1)*t(k+2)+(2*R*k+2*R)*t(k+1)+(R^2*(k-1)+2*R^2)*t(k)-(2*R*k*(k+1)+2*R*(k+1))*y(k+2)+(-R^2*(k-1)*k+D0*S-b*S*n-2*R^2*k)*y(k+1)+(D1*S-2*b*S*n*R)*y(k)+(D2*S-b*S*n*R^2)*y(k-1))/((1+DgS)*(k+1)*(k+2));
        t(k)=(-(4*S*R*k*(k+1)+4*S*R*(k+1))*t(k+2)+(1+FS-b*S*ra-6*S*R^2*(k-1)*k-12*S*R^2*k)*t(k+1)+(-4*S*R^3*(k-2)*(k-1)+2*R-4*b*S*ra*R-12*S*R^3*(k-1))*t(k)+(-S*R^4*(k-3)*(k-2)+R^2-6*b*S*ra*R^2-4*S*R^4*(k-2))*t(k-1)-(4*b*S*ra*R^3)*t(k-2)-(b*S*ra*R^4)*t(k-3)-(1+FS)*(k+1)*y(k+2)-2*R*k*y(k+1)-R^2*(k-1)*y(k))/(S*(k+1)*(k+2));
end 

Now consider execution for the first iteration, when k = 3.

y(t) uses: t(k+2), t(k+1), t(k), y(k+2), y(k+1), y(k), y(k-1) . The smallest index used there is y(k-1) which would be y(3-1) --> y(2) which would be fine. It is a bit questionable about whether it should be using y(k) to define y(k) but doing so would result in any subscript problems as y was pre-initialized.

t(k) uses: t(k+2), t(k+1), t(k), t(k-1), t(k-2), t(k-3), y(k+2), y(k+1), y(k) . The smallest index used there is t(k-3) which would be y(3-3) --> y(0) which is a problem.

It is not possible to define t(0) in MATLAB.

If your posted code is different than your actual code and your actual code does not use t(k-3) or your actual code starts with k = 4, then we need to see your actual code to debug it. But for your posted code, you have the problem that t(k-3) is t(0) when k = 3.

[Andres Cordoba Ordoñez]

function [fi, fi1, Ji, Ji1] = funcion(L,r1,s1,t1,C0,C1,C2,C3,FS,b,D0,D1,D2,AstG,S,R,DgS,ra,n)
% syms x  ; % Definimos las variables simbólicas
%definimos las condiciones iniciales
m=100;
% mi=m+3; %# de terminos
% % Definimos una variable "y" como un vector de longitud m n función de la variable simbólica L.
y = sym(zeros(1,m));
t = sym(zeros(1,m));
    y(1)=C0;
    y(2)=C1;
    t(1)=C2;
    t(2)=C3;
    y(3)=((1+FS)*t(2)+(2*R)*t(1)-(2*R*(1))*y(2)+(D0*S-b*S*n)*y(1)+L*r1/AstG)/((1+DgS)*2);
    t(3)=(-(4*S*R*(1))*t(2)+(1+FS-b*S*ra)*t(1)-(1+FS)*(1)*y(2))/(S*2);
    y(4)=((1+FS)*(2)*t(3)+(2*R*1+2*R)*t(2)+(2*R^2)*t(1)-(2*R*1*(2)+2*R*(2))*y(3)+(D0*S-b*S*n-2*R^2*1)*y(2)+(D1*S-2*b*S*n*R)*y(1)+L^2*s1/AstG)/((1+DgS)*6);
    t(4)=(-(4*S*R*1*(2)+4*S*R*(2))*t(3)+(1+FS-b*S*ra-12*S*R^2*1)*t(2)+(2*R-4*b*S*ra*R)*t(1)-(1+FS)*(2)*y(3)-2*R*1*y(2))/(S*6);
    y(5)=((1+FS)*(3)*t(4)+(2*R*2+2*R)*t(3)+(R^2*(1)+2*R^2)*t(2)-(2*R*2*(3)+2*R*(3))*y(4)+(-R^2*(1)*2+D0*S-b*S*n-2*R^2*2)*y(3)+(D1*S-2*b*S*n*R)*y(2)+(D2*S-b*S*n*R^2)*y(1)+(t1*L^3)/AstG)/((1+DgS)*12);
    t(5)=(-(4*S*R*2*(3)+4*S*R*(3))*t(4)+(1+FS-b*S*ra-6*S*R^2*(1)*2-12*S*R^2*2)*t(3)+(2*R-4*b*S*ra*R-12*S*R^3*(1))*t(2)+(R^2-6*b*S*ra*R^2)*t(1)-(1+FS)*(3)*y(4)-2*R*2*y(3)-R^2*(1)*y(2))/(S*12);
    y(6)=((1+FS)*(4)*t(5)+(2*R*3+2*R)*t(4)+(R^2*(2)+2*R^2)*t(3)-(2*R*3*(4)+2*R*(4))*y(5)+(-R^2*(2)*3+D0*S-b*S*n-2*R^2*3)*y(4)+(D1*S-2*b*S*n*R)*y(3)+(D2*S-b*S*n*R^2)*y(2))/((1+DgS)*20);
    t(6)=(-(4*S*R*3*(4)+4*S*R*(4))*t(5)+(1+FS-b*S*ra-6*S*R^2*(2)*3-12*S*R^2*3)*t(4)+(-4*S*R^3*(1)*(2)+2*R-4*b*S*ra*R-12*S*R^3*(2))*t(3)+(R^2-6*b*S*ra*R^2-4*S*R^4*(1))*t(2)-(4*b*S*ra*R^3)*t(1)-(1+FS)*(4)*y(5)-2*R*3*y(4)-R^2*(2)*y(3))/(S*20);
    
for k = 4:m
    
    y(k+3)=((1+FS)*(k+1)*t(k+2)+(2*R*k+2*R)*t(k+1)+(R^2*(k-1)+2*R^2)*t(k)-(2*R*k*(k+1)+2*R*(k+1))*y(k+2)+(-R^2*(k-1)*k+D0*S-b*S*n-2*R^2*k)*y(k+1)+(D1*S-2*b*S*n*R)*y(k)+(D2*S-b*S*n*R^2)*y(k-1))/((1+DgS)*(k+1)*(k+2));
    
    t(k+3)=(-(4*S*R*k*(k+1)+4*S*R*(k+1))*t(k+2)+(1+FS-b*S*ra-6*S*R^2*(k-1)*k-12*S*R^2*k)*t(k+1)+(-4*S*R^3*(k-2)*(k-1)+2*R-4*b*S*ra*R-12*S*R^3*(k-1))*t(k)+(-S*R^4*(k-3)*(k-2)+R^2-6*b*S*ra*R^2-4*S*R^4*(k-2))*t(k-1)-(4*b*S*ra*R^3)*t(k-2)-(b*S*ra*R^4)*t(k-3)-(1+FS)*(k+1)*y(k+2)-2*R*k*y(k+1)-R^2*(k-1)*y(k))/(S*(k+1)*(k+2));
end 
j = 0:m+2;
fi = sum(y);
fi1 = y * j';
Ji = sum(t);
Ji1 = t * j';
                                  
end     

[Walter Roberson]

I am not getting an error? (Note: I had to reduce m in order to finish execution within the Answers time limit.)

syms(str2sym('[L,r1,s1,t1,C0,C1,C2,C3,FS,b,D0,D1,D2,AstG,S,R,DgS,ra,n]'))
funcion(L,r1,s1,t1,C0,C1,C2,C3,FS,b,D0,D1,D2,AstG,S,R,DgS,ra,n)

ans = 

function [fi, fi1, Ji, Ji1] = funcion(L,r1,s1,t1,C0,C1,C2,C3,FS,b,D0,D1,D2,AstG,S,R,DgS,ra,n)
% syms x  ; % Definimos las variables simbólicas
%definimos las condiciones iniciales
m = 6; %100;
% mi=m+3; %# de terminos
% % Definimos una variable "y" como un vector de longitud m n función de la variable simbólica L.
y = sym(zeros(1,m));
t = sym(zeros(1,m));
    y(1)=C0;
    y(2)=C1;
    t(1)=C2;
    t(2)=C3;
    y(3)=((1+FS)*t(2)+(2*R)*t(1)-(2*R*(1))*y(2)+(D0*S-b*S*n)*y(1)+L*r1/AstG)/((1+DgS)*2);
    t(3)=(-(4*S*R*(1))*t(2)+(1+FS-b*S*ra)*t(1)-(1+FS)*(1)*y(2))/(S*2);
    y(4)=((1+FS)*(2)*t(3)+(2*R*1+2*R)*t(2)+(2*R^2)*t(1)-(2*R*1*(2)+2*R*(2))*y(3)+(D0*S-b*S*n-2*R^2*1)*y(2)+(D1*S-2*b*S*n*R)*y(1)+L^2*s1/AstG)/((1+DgS)*6);
    t(4)=(-(4*S*R*1*(2)+4*S*R*(2))*t(3)+(1+FS-b*S*ra-12*S*R^2*1)*t(2)+(2*R-4*b*S*ra*R)*t(1)-(1+FS)*(2)*y(3)-2*R*1*y(2))/(S*6);
    y(5)=((1+FS)*(3)*t(4)+(2*R*2+2*R)*t(3)+(R^2*(1)+2*R^2)*t(2)-(2*R*2*(3)+2*R*(3))*y(4)+(-R^2*(1)*2+D0*S-b*S*n-2*R^2*2)*y(3)+(D1*S-2*b*S*n*R)*y(2)+(D2*S-b*S*n*R^2)*y(1)+(t1*L^3)/AstG)/((1+DgS)*12);
    t(5)=(-(4*S*R*2*(3)+4*S*R*(3))*t(4)+(1+FS-b*S*ra-6*S*R^2*(1)*2-12*S*R^2*2)*t(3)+(2*R-4*b*S*ra*R-12*S*R^3*(1))*t(2)+(R^2-6*b*S*ra*R^2)*t(1)-(1+FS)*(3)*y(4)-2*R*2*y(3)-R^2*(1)*y(2))/(S*12);
    y(6)=((1+FS)*(4)*t(5)+(2*R*3+2*R)*t(4)+(R^2*(2)+2*R^2)*t(3)-(2*R*3*(4)+2*R*(4))*y(5)+(-R^2*(2)*3+D0*S-b*S*n-2*R^2*3)*y(4)+(D1*S-2*b*S*n*R)*y(3)+(D2*S-b*S*n*R^2)*y(2))/((1+DgS)*20);
    t(6)=(-(4*S*R*3*(4)+4*S*R*(4))*t(5)+(1+FS-b*S*ra-6*S*R^2*(2)*3-12*S*R^2*3)*t(4)+(-4*S*R^3*(1)*(2)+2*R-4*b*S*ra*R-12*S*R^3*(2))*t(3)+(R^2-6*b*S*ra*R^2-4*S*R^4*(1))*t(2)-(4*b*S*ra*R^3)*t(1)-(1+FS)*(4)*y(5)-2*R*3*y(4)-R^2*(2)*y(3))/(S*20);
    
for k = 4:m
    
    y(k+3)=((1+FS)*(k+1)*t(k+2)+(2*R*k+2*R)*t(k+1)+(R^2*(k-1)+2*R^2)*t(k)-(2*R*k*(k+1)+2*R*(k+1))*y(k+2)+(-R^2*(k-1)*k+D0*S-b*S*n-2*R^2*k)*y(k+1)+(D1*S-2*b*S*n*R)*y(k)+(D2*S-b*S*n*R^2)*y(k-1))/((1+DgS)*(k+1)*(k+2));
    
    t(k+3)=(-(4*S*R*k*(k+1)+4*S*R*(k+1))*t(k+2)+(1+FS-b*S*ra-6*S*R^2*(k-1)*k-12*S*R^2*k)*t(k+1)+(-4*S*R^3*(k-2)*(k-1)+2*R-4*b*S*ra*R-12*S*R^3*(k-1))*t(k)+(-S*R^4*(k-3)*(k-2)+R^2-6*b*S*ra*R^2-4*S*R^4*(k-2))*t(k-1)-(4*b*S*ra*R^3)*t(k-2)-(b*S*ra*R^4)*t(k-3)-(1+FS)*(k+1)*y(k+2)-2*R*k*y(k+1)-R^2*(k-1)*y(k))/(S*(k+1)*(k+2));
end 
j = 0:m+2;
fi = sum(y);
fi1 = y * j';
Ji = sum(t);
Ji1 = t * j';
                                  
end     

[Andres Cordoba Ordoñez]

the question is; how can improve my code in order to be faster, if a tried to change the solution from syms to double, it lost precision and the answer is wrong.

[Walter Roberson]

According to my tests, a pure numeric solution (no syms at all) is pretty much 100 times faster, and is accurate to one part in

Side note: my tests with lower m values showed that at least for these random inputs, the calculations are divergent, increasing in value moderately rapidly with larger m values.

To do further timing and accuracy tests, we will need your inputs.

format long g
rng(655321)
values = num2cell(randn(1, 19) * 10);
[L,r1,s1,t1,C0,C1,C2,C3,FS,b,D0,D1,D2,AstG,S,R,DgS,ra,n] = values{:};

tic
[OUT1, OUT2, OUT3, OUT4] = funcion(L,r1,s1,t1,C0,C1,C2,C3,FS,b,D0,D1,D2,AstG,S,R,DgS,ra,n);
OUT1n = double(OUT1)

OUT1n = 
     -7.24193443872753e+17

OUT2n = double(OUT2)

OUT2n = 
     -7.41640401604103e+19

OUT3n = double(OUT3)

OUT3n = 
       3.1885102761081e+20

OUT4n = double(OUT4)

OUT4n = 
      3.26527272151145e+22

toc

Elapsed time is 4.381793 seconds.

tic
[OUT1_num, OUT2_num, OUT3_num, OUT4_num] = funcion_num(L,r1,s1,t1,C0,C1,C2,C3,FS,b,D0,D1,D2,AstG,S,R,DgS,ra,n);
OUT1n_num = double(OUT1_num)

OUT1n_num = 
     -7.24193443871126e+17

OUT2n_num = double(OUT2_num)

OUT2n_num = 
     -7.41640401602437e+19

OUT3n_num = double(OUT3_num)

OUT3n_num = 
       3.1885102761008e+20

OUT4n_num = double(OUT4_num)

OUT4n_num = 
      3.26527272150397e+22

toc

Elapsed time is 0.041566 seconds.

OUT1n - OUT1n_num, ans/OUT1n

ans = 
    -1626496

ans = 
      2.24594134862921e-12

OUT2n - OUT2n_num, ans/OUT2n

ans = 
  -166559744

ans = 
      2.24582889011637e-12

OUT3n - OUT3n_num, ans/OUT3n

ans = 
   730136576

ans = 
      2.28989877019059e-12

OUT4n - OUT4n_num, ans/OUT4n

ans = 
               74767663104

ans = 
      2.28978310483638e-12

values_sym = num2cell(sym(values));
[L,r1,s1,t1,C0,C1,C2,C3,FS,b,D0,D1,D2,AstG,S,R,DgS,ra,n] = values{:};

tic
[OUT1s, OUT2s, OUT3s, OUT4s] = funcion(L,r1,s1,t1,C0,C1,C2,C3,FS,b,D0,D1,D2,AstG,S,R,DgS,ra,n);
OUT1ns = double(OUT1s)

OUT1ns = 
     -7.24193443872753e+17

OUT2ns = double(OUT2s)

OUT2ns = 
     -7.41640401604103e+19

OUT3ns = double(OUT3s)

OUT3ns = 
       3.1885102761081e+20

OUT4ns = double(OUT4s)

OUT4ns = 
      3.26527272151145e+22

toc

Elapsed time is 3.313651 seconds.

OUT1n - OUT1ns

ans = 
     0

OUT2n - OUT2ns

ans = 
     0

OUT3n - OUT3ns

ans = 
     0

OUT4n - OUT4ns

ans = 
     0

function [fi, fi1, Ji, Ji1] = funcion(L,r1,s1,t1,C0,C1,C2,C3,FS,b,D0,D1,D2,AstG,S,R,DgS,ra,n)
% syms x  ; % Definimos las variables simbólicas
%definimos las condiciones iniciales
m = 100;
% mi=m+3; %# de terminos
% % Definimos una variable "y" como un vector de longitud m n función de la variable simbólica L.
y = sym(zeros(1,m));
t = sym(zeros(1,m));
    y(1)=C0;
    y(2)=C1;
    t(1)=C2;
    t(2)=C3;
    y(3)=((1+FS)*t(2)+(2*R)*t(1)-(2*R*(1))*y(2)+(D0*S-b*S*n)*y(1)+L*r1/AstG)/((1+DgS)*2);
    t(3)=(-(4*S*R*(1))*t(2)+(1+FS-b*S*ra)*t(1)-(1+FS)*(1)*y(2))/(S*2);
    y(4)=((1+FS)*(2)*t(3)+(2*R*1+2*R)*t(2)+(2*R^2)*t(1)-(2*R*1*(2)+2*R*(2))*y(3)+(D0*S-b*S*n-2*R^2*1)*y(2)+(D1*S-2*b*S*n*R)*y(1)+L^2*s1/AstG)/((1+DgS)*6);
    t(4)=(-(4*S*R*1*(2)+4*S*R*(2))*t(3)+(1+FS-b*S*ra-12*S*R^2*1)*t(2)+(2*R-4*b*S*ra*R)*t(1)-(1+FS)*(2)*y(3)-2*R*1*y(2))/(S*6);
    y(5)=((1+FS)*(3)*t(4)+(2*R*2+2*R)*t(3)+(R^2*(1)+2*R^2)*t(2)-(2*R*2*(3)+2*R*(3))*y(4)+(-R^2*(1)*2+D0*S-b*S*n-2*R^2*2)*y(3)+(D1*S-2*b*S*n*R)*y(2)+(D2*S-b*S*n*R^2)*y(1)+(t1*L^3)/AstG)/((1+DgS)*12);
    t(5)=(-(4*S*R*2*(3)+4*S*R*(3))*t(4)+(1+FS-b*S*ra-6*S*R^2*(1)*2-12*S*R^2*2)*t(3)+(2*R-4*b*S*ra*R-12*S*R^3*(1))*t(2)+(R^2-6*b*S*ra*R^2)*t(1)-(1+FS)*(3)*y(4)-2*R*2*y(3)-R^2*(1)*y(2))/(S*12);
    y(6)=((1+FS)*(4)*t(5)+(2*R*3+2*R)*t(4)+(R^2*(2)+2*R^2)*t(3)-(2*R*3*(4)+2*R*(4))*y(5)+(-R^2*(2)*3+D0*S-b*S*n-2*R^2*3)*y(4)+(D1*S-2*b*S*n*R)*y(3)+(D2*S-b*S*n*R^2)*y(2))/((1+DgS)*20);
    t(6)=(-(4*S*R*3*(4)+4*S*R*(4))*t(5)+(1+FS-b*S*ra-6*S*R^2*(2)*3-12*S*R^2*3)*t(4)+(-4*S*R^3*(1)*(2)+2*R-4*b*S*ra*R-12*S*R^3*(2))*t(3)+(R^2-6*b*S*ra*R^2-4*S*R^4*(1))*t(2)-(4*b*S*ra*R^3)*t(1)-(1+FS)*(4)*y(5)-2*R*3*y(4)-R^2*(2)*y(3))/(S*20);
    
for k = 4:m
    
    y(k+3)=((1+FS)*(k+1)*t(k+2)+(2*R*k+2*R)*t(k+1)+(R^2*(k-1)+2*R^2)*t(k)-(2*R*k*(k+1)+2*R*(k+1))*y(k+2)+(-R^2*(k-1)*k+D0*S-b*S*n-2*R^2*k)*y(k+1)+(D1*S-2*b*S*n*R)*y(k)+(D2*S-b*S*n*R^2)*y(k-1))/((1+DgS)*(k+1)*(k+2));
    
    t(k+3)=(-(4*S*R*k*(k+1)+4*S*R*(k+1))*t(k+2)+(1+FS-b*S*ra-6*S*R^2*(k-1)*k-12*S*R^2*k)*t(k+1)+(-4*S*R^3*(k-2)*(k-1)+2*R-4*b*S*ra*R-12*S*R^3*(k-1))*t(k)+(-S*R^4*(k-3)*(k-2)+R^2-6*b*S*ra*R^2-4*S*R^4*(k-2))*t(k-1)-(4*b*S*ra*R^3)*t(k-2)-(b*S*ra*R^4)*t(k-3)-(1+FS)*(k+1)*y(k+2)-2*R*k*y(k+1)-R^2*(k-1)*y(k))/(S*(k+1)*(k+2));
end 
j = 0:m+2;
fi = sum(y);
fi1 = y * j';
Ji = sum(t);
Ji1 = t * j';
                                  
end 

function [fi, fi1, Ji, Ji1] = funcion_num(L,r1,s1,t1,C0,C1,C2,C3,FS,b,D0,D1,D2,AstG,S,R,DgS,ra,n)
% syms x  ; % Definimos las variables simbólicas
%definimos las condiciones iniciales
m = 100;
% mi=m+3; %# de terminos
% % Definimos una variable "y" como un vector de longitud m n función de la variable simbólica L.
y = (zeros(1,m));
t = (zeros(1,m));
    y(1)=C0;
    y(2)=C1;
    t(1)=C2;
    t(2)=C3;
    y(3)=((1+FS)*t(2)+(2*R)*t(1)-(2*R*(1))*y(2)+(D0*S-b*S*n)*y(1)+L*r1/AstG)/((1+DgS)*2);
    t(3)=(-(4*S*R*(1))*t(2)+(1+FS-b*S*ra)*t(1)-(1+FS)*(1)*y(2))/(S*2);
    y(4)=((1+FS)*(2)*t(3)+(2*R*1+2*R)*t(2)+(2*R^2)*t(1)-(2*R*1*(2)+2*R*(2))*y(3)+(D0*S-b*S*n-2*R^2*1)*y(2)+(D1*S-2*b*S*n*R)*y(1)+L^2*s1/AstG)/((1+DgS)*6);
    t(4)=(-(4*S*R*1*(2)+4*S*R*(2))*t(3)+(1+FS-b*S*ra-12*S*R^2*1)*t(2)+(2*R-4*b*S*ra*R)*t(1)-(1+FS)*(2)*y(3)-2*R*1*y(2))/(S*6);
    y(5)=((1+FS)*(3)*t(4)+(2*R*2+2*R)*t(3)+(R^2*(1)+2*R^2)*t(2)-(2*R*2*(3)+2*R*(3))*y(4)+(-R^2*(1)*2+D0*S-b*S*n-2*R^2*2)*y(3)+(D1*S-2*b*S*n*R)*y(2)+(D2*S-b*S*n*R^2)*y(1)+(t1*L^3)/AstG)/((1+DgS)*12);
    t(5)=(-(4*S*R*2*(3)+4*S*R*(3))*t(4)+(1+FS-b*S*ra-6*S*R^2*(1)*2-12*S*R^2*2)*t(3)+(2*R-4*b*S*ra*R-12*S*R^3*(1))*t(2)+(R^2-6*b*S*ra*R^2)*t(1)-(1+FS)*(3)*y(4)-2*R*2*y(3)-R^2*(1)*y(2))/(S*12);
    y(6)=((1+FS)*(4)*t(5)+(2*R*3+2*R)*t(4)+(R^2*(2)+2*R^2)*t(3)-(2*R*3*(4)+2*R*(4))*y(5)+(-R^2*(2)*3+D0*S-b*S*n-2*R^2*3)*y(4)+(D1*S-2*b*S*n*R)*y(3)+(D2*S-b*S*n*R^2)*y(2))/((1+DgS)*20);
    t(6)=(-(4*S*R*3*(4)+4*S*R*(4))*t(5)+(1+FS-b*S*ra-6*S*R^2*(2)*3-12*S*R^2*3)*t(4)+(-4*S*R^3*(1)*(2)+2*R-4*b*S*ra*R-12*S*R^3*(2))*t(3)+(R^2-6*b*S*ra*R^2-4*S*R^4*(1))*t(2)-(4*b*S*ra*R^3)*t(1)-(1+FS)*(4)*y(5)-2*R*3*y(4)-R^2*(2)*y(3))/(S*20);
    
for k = 4:m
    
    y(k+3)=((1+FS)*(k+1)*t(k+2)+(2*R*k+2*R)*t(k+1)+(R^2*(k-1)+2*R^2)*t(k)-(2*R*k*(k+1)+2*R*(k+1))*y(k+2)+(-R^2*(k-1)*k+D0*S-b*S*n-2*R^2*k)*y(k+1)+(D1*S-2*b*S*n*R)*y(k)+(D2*S-b*S*n*R^2)*y(k-1))/((1+DgS)*(k+1)*(k+2));
    
    t(k+3)=(-(4*S*R*k*(k+1)+4*S*R*(k+1))*t(k+2)+(1+FS-b*S*ra-6*S*R^2*(k-1)*k-12*S*R^2*k)*t(k+1)+(-4*S*R^3*(k-2)*(k-1)+2*R-4*b*S*ra*R-12*S*R^3*(k-1))*t(k)+(-S*R^4*(k-3)*(k-2)+R^2-6*b*S*ra*R^2-4*S*R^4*(k-2))*t(k-1)-(4*b*S*ra*R^3)*t(k-2)-(b*S*ra*R^4)*t(k-3)-(1+FS)*(k+1)*y(k+2)-2*R*k*y(k+1)-R^2*(k-1)*y(k))/(S*(k+1)*(k+2));
end 
j = 0:m+2;
fi = sum(y);
fi1 = y * j';
Ji = sum(t);
Ji1 = t * j';
                                  
end