Fractional order algebraic equations

Hi, I am having a system of two fractional differential equations and an algebraic equation. I would like to solve that system and plot the time series. We defining the algebraic equation in matlab, an error shown as an invalid matlab expression. The code is as follows function [T, Y]=FOAE(parameters, orders, TSim, Y0) %

% time step:
h=0.01; 
% number of calculated mesh points:
n=round(TSim/h);
%orders of derivatives, respectively:
q1=orders(1); q2=orders(2); q3=orders(3); 
% constants of Duffing's system:
 r=parameters(1); k=parameters(2); beta=parameters(3);sigma=parameters(4); a=parameters(5); p=parameters(6); c=parameters(7); m=parameters(8);
% binomial coefficients calculation:
cp1=1; cp2=1; cp3=1;
for j=1:n
    c1(j)=(1-(1+q1)/j)*cp1;
    c2(j)=(1-(1+q2)/j)*cp2;
    c3(j)=(1-(1+q3)/j)*cp3;
    cp1=c1(j); cp2=c2(j); cp3=c3(j);    
end
% initial conditions setting:
x(1)=Y0(1); y(1)=Y0(2); z(1)=Y0(3);
% calculation of phase portraits /numerical solution/:
for i=2:n
    x(i)=(r*x(i-1)*(1-x(i-1)/k)-beta*x(i-1)*y(i-1)/(1+sigma*x(i-1)))*h^q1 - memo(x, c1, i);
    y(i)=(beta*x(i-1)*y(i-1)/(1+sigma*x(i-1))-a*y(i-1)-y(i-1)*z(i-1))*h^q2 - memo(y, c2, i);
    0   =(z(i-1)*(p*y(i-1)-c)-m)*h^q3 - memo(z, c3, i);
end
for j=1:n
    Y(j,1)=x(j);
    Y(j,2)=y(j);
    Y(j,3)=z(j);
end
T=0:h:TSim;
%%%%%%%%%%%%%%%%%
[t, y]=FOAE([0.2 5 0.2 0.001 0.2 1.5 1 -0.0001], [  0.8    0.8  0 ], 100, [1.3 0.4 0.00025]);
 figure(1)
plot(y(:,1), y(:,2), 'k','Linewidth',3);
fsize=15;
xlabel('x_1(t)','FontSize',fsize); ylabel('x_2(t)','FontSize',fsize); grid;

%%%%%%%%%%%%%%%%%%%%%% Any help is really appreciated. Thank you very much in advance. Sanaa Moussa

[Matt J]

Please indent your lines of code to make them split across lines readably.

[Walter Roberson]

 0 =(z(i-1)*(p*y(i-1)-c)-m)*h^q3 - memo(z, c3, i);

Is not a valid line. It attempts to assign an expression to the constant 0