Solution :
cs=zeros(9001);
ca=zeros(9001);
cp=zeros(9001);
psi=zeros(9001);
t=[0:0.1:900];
cs(1)=0.5;
ce(1)=0.001;
cp(1)=0;
ca(1)=0;
psi(1)=0;
for i=1:1:9000
cs(i+1)=cs(i)-0.1*((0.015*cs(i))/(5.53+cs(i)));
cp(i+1)=cp(i)+0.1*((0.015*cs(i))/(5.53+cs(i))-0.0026*cp(i));
ca(i+1)=ca(i)+0.1*0.0026*cp(i);
psi(i+1)=((cp(i+1)-cp(i)))/((cs(i)-cs(i+1)));
end
plot(t,cs,t,cp,t,ca);
plot(t,psi);