this is the program
plsolve
t=1:20;
a(1)=1;
b(1)=0;
c(1)=1;
sigma=10;
r=28;
d=(8/3);
sumA=0;
sumC=0;
sumB=0;
sumS=0;
sumZ=0;
for i=1:30
for p=1:29;
sumS=sumS+(a(p).*c(i-p));
sumZ=sumZ+(a(p).*b(i-p));
a(i+1)=((sigma)/(i+1))*(a(i)-b(i));
b(i+1)=(1/(i+1))*(r*a(i)-b(i)-sumS);
c(i+1)=(1/(i+1))*(sumZ-d*c(i));
sumA=sumA+a(i).*(t^i);
sumB=sumB+b(i).*(t^i);
sumC=sumC+c(i).*(t^1);
end
end
plot(t,sumB)
