Vedat Erturk I'll try to help: give you some ideas to try. Using an array for: α1=1;α2=1;α3=1;α4=1;α5=1; α={1,1,1,1,1} and referencing as: α[[I]] could help. Maybe better as a[i_].Other variables could also be assigned to an array/ given a number instead of a variable name. In most model Intel based machines :ParallelTable[] and ParallelSum[] speed up calculations by the number of processes your machine allows. Integrations might be better as NIntegrate for time saving. Do Loops have ever been slow in Mathematica. Replacing the loops with Block[] might help as sometimes used in fractal programs. Compiling the Do loop section might also speed up the calculation process. I know the Parallel idea works for at least 10 places in your program. In all: buy yourself the fastest computer you can get (with as many processors as you can) with the most built in memory (16 gb or more): then the slowness of Mathematica will fade, LOL. Good luck in your real world programming. Roger Lee Bagula

This is the current state of my modifications

ClearAll["Global`*"];
(*modified code to eliminate subscripts and underscores*)
(*modified code changing z[] to z2f[] and z3f[]*)
α1=1;α2=1;α3=1;α4=1;α5=1;
g0=0.9;b10=0.03;b11=0.001;Topt=24;gama0=150;ds0=4;A0=0.2;
w=2.10;z1=3.5;Oc=1.10;d1=0.66;d=0.4;gama1=2;d2=0.2;z=0.0019;
β=1.024;γ=4;z0=10.10;T10=14.50;delta1=0.1;T0=20;
b=1.30;Tmaks=35;K=2;M=2;T=60.05;
x[t_,0,0]:=10;h1[t_,0,0]:=28;h2[t_,0,0]:=1;y[t_,0,0]:=1.2;z3f[t_,0,0]:=0.25;
x[ λ_,m_]:=Sum[ x[t,m,i]*λ^i,{i,0,K}]; 
h1[λ_,m_]:=Sum[h1[t,m,i]*λ^i,{i,0,K}]; 
h2[λ_,m_]:=Sum[h2[t,m,i]*λ^i,{i,0,K}]; 
y[ λ_,m_]:=Sum[ y[t,m,i]*λ^i,{i,0,K}];
z2f[λ_,m_]:=Sum[z3f[t,m,i]*λ^i,{i,0,K}];
N1[λ_,m_]:=g0*(Exp[-b*(h1[λ,m]-Topt)/(Tmaks-Topt)]+(z2f[λ,m]-(b10+b11*
(h1[λ,m]-Topt)))/(z2f[λ,m]+1))*x[λ,m]-g0*x[λ,m]^2/gama0;
N2[λ_,m_]:=w*h2[λ,m]-z1*(h1[λ,m]-T10)+γ*(z0-y[λ,m]);
N3[λ_,m_]:=A0-delta1*h2[λ,m];
N4[λ_,m_]:=Oc-d1*y[λ,m]-d*y[λ,m]*h2[λ,m];
N5[λ_,m_]:=gama1*β^(h1[λ,m]-T0)*(ds0/(1+h1[λ,m]-Topt)-z2f[λ,m])-
d2*z2f[λ,m]*x[λ,m]-z*(h1[λ,m]-Topt);
A1[m_,0]:=N1[0,m];A2[m_,0]:=N2[0,m];A3[m_,0]:=N3[0,m];A4[m_,0]:=N4[0,m];A5[m_,0]:=N5[0,m];
Do[
  For[i=0,i<K,i++,
    A1[m,i]=Simplify[1/i!*D[N1[λ,m],{λ,i}]]/.λ->0;
    Print["{i,m,derivative}=",{i,m,A1[m,i]}];
    A2[m,i]=Simplify[1/i!*D[N2[λ,m],{λ,i}]]/.λ->0;
    A3[m,i]=Simplify[1/i!*D[N3[λ,m],{λ,i}]]/.λ->0;
    A4[m,i]=Simplify[1/i!*D[N4[λ,m],{λ,i}]]/.λ->0;
    A5[m,i]=Simplify[1/i!*D[N5[λ,m],{λ,i}]]/.λ->0
  ];
  Do[
    x[ t_,m,k]=1/Gamma[α1]*Integrate[(t-τ)^(α1-1)*A1[m,k-1],{τ,0,t}];
    Print["{k,m,integral}=",{k,m,Simplify[1/Gamma[\[Alpha]1]*Integrate[(t-\[Tau])^(\[Alpha]1-1)*A1[m,k-1],{\[Tau],0,t}]]}];
    h1[t_,m,k]=1/Gamma[α2]*Integrate[(t-τ)^(α2-1)*A2[m,k-1],{τ,0,t}];
    h2[t_,m,k]=1/Gamma[α3]*Integrate[(t-τ)^(α3-1)*A3[m,k-1],{τ,0,t}];
    y[ t_,m,k]=1/Gamma[α4]*Integrate[(t-τ)^(α4-1)*A4[m,k-1],{τ,0,t}];
    z3f[t_,m,k]=1/Gamma[α5]*Integrate[(t-τ)^(α5-1)*A5[m,k-1],{τ,0,t}],
  {k,1,K}];
  f1[t_,m_]:=Sum[x[ t-(T/M)*m,m,i],{i,0,K}];
  f2[t_,m_]:=Sum[h1[t-(T/M)*m,m,i],{i,0,K}];
  f3[t_,m_]:=Sum[h2[t-(T/M)*m,m,i],{i,0,K}];
  f4[t_,m_]:=Sum[y[ t-(T/M)*m,m,i],{i,0,K}];
  f5[t_,m_]:=Sum[z3f[t-(T/M)*m,m,i],{i,0,K}];
  x[ t_,m+1,0]=f1[(T/M)*(m+1),m];
  h1[t_,m+1,0]=f2[(T/M)*(m+1),m];
  h2[t_,m+1,0]=f3[(T/M)*(m+1),m];
  y[ t_,m+1,0]=f4[(T/M)*(m+1),m];
  z3f[t_,m+1,0]=f5[(T/M)*(m+1),m],
{m,0,M-1}];
(*X[ t_]=X[ t]=Piecewise[Table[{f1[t,m],((T/M)*(m))<=t<((T/M)*(m+1))},{m,0,M-1}]];
  H1[t_]=H1[t]=Piecewise[Table[{f2[t,m],((T/M)*(m))<=t<((T/M)*(m+1))},{m,0,M-1}]];
  H2[t_]=H2[t]=Piecewise[Table[{f3[t,m],((T/M)*(m))<=t<((T/M)*(m+1))},{m,0,M-1}]];
  Y[ t_]=Y[ t]=Piecewise[Table[{f4[t,m],((T/M)*(m))<=t<((T/M)*(m+1))},{m,0,M-1}]];
  Z[ t_]=Z[ t]=Piecewise[Table[{f5[t,m],((T/M)*(m))<=t<((T/M)*(m+1))},{m,0,M-1}]];
  s1=Table[{t,X[t]},{t,0,T}]
  s1MSADM=ListLinePlot[s1,PlotRange->All,PlotStyle->Black,AxesLabel->{"t","N(t)"}]*)

And the output of that is

{i,m,derivative}={0,0,6.5649}
{i,m,derivative}={1,0,0.59649 x[t,0,1]+9.(-0.0736627 h1[t,0,1]-0.13824([0.+z3f[t,0,1])+0.8(0.-0.001 h1[t,0,1]+z3f[t,0,1]))}
{i,m,integral}={2,0,6.56490*t}
{i,m,integral}={2,0,0+15.2130*t^2}
{i,m,derivative}={0,1,0}
{i,m,derivative}={1,1,0}
{i,m,integral}={2,1,0}
{i,m,integral}={2,1,0}
{i,m,derivative}={0,2,0}
{i,m,derivative}={1,2,0}
{i,m,integral}={2,2,0}
{i,m,integral}={2,2,0}

I believe what that is showing is that you might not need all the machinery and horsepower to be doing derivatives and integrals to construct your Piecewise function. If that is true then this may result in substantial increases in the speed of your code, which was your stated goal. If you look at the code with the idea of finding further simplifications in the way it calculates the final result then I think it might be possible to simplify and speed this up even further.

This still needs more work and testing and verification. Please check this carefully to try to make certain that I have not introduced any errors.

Note: There are tiny changes between the code you first posted and the notebook you later posted. I believe my code is based on the original post. I also think I have probably made mistakes translating your code to mine and that my calculations probably do not match yours. If you can point out specific differences for i,m,k and the derivative or integral, what you get and I should get then I will try to track down my errors and see if I can fix those. It is also likely that the diagnostic prints that I inserted could be improved to provide better information. Please make any suggestions.

Thanks

I should be able to get the K at least 10.

Hello Bill, Thank you for your reply.only x is double subscripted. K is 5,10,20,30,etc. that is,it is not so big.some times 10 is even enough. Now, T/M is just 0.1. I should be able to choose T/M as small as I want. For example T/M=0.01,0.001,etc. Hello again Bill, I hope this helps you.