V=5.07*10^10;
Q=2.085*10^9;
Dt=0.5;  %Dt
R=10950
T=V/Q;     %XRONOS PARAMONIS
L=201*10^9; %ug/day 
W1=0.9*L;
W2=0.1*L;
N=zeros(1,R+1);
P=zeros(1,R+1);
No=zeros(1,R+1);
N(1,1)=10.215;    %ug/L
P(1,1)=2.73;     %ug/L
No(1,1)=1.135;    %ug/L
Kg=0.04;         %L/ug/day
Ko=0.9;      %1/day
K1=0.4067;       %1/day
Ks=0.001;       %1/day
for i=1:R
    DN(1,i)=W1/V-N(1,i)/T+Ko*No(1,i)-Kg*N(1,i)*P(1,i);
    DP(1,i)=Kg*N(1,i)*P(1,i)-(K1+Ks)*P(1,i)-P(1,i)/T;
    DNo(1,i)=W2/V-No(1,i)/T+K1*P(1,i)-(Ks+Ko)*No(1,i);
    N(1,i+1)=N(1,i)+DN(1,i)*Dt;
    P(1,i+1)=P(1,i)+DP(1,i)*Dt;
    No(1,i+1)=No(1,i)+DNo(1,i)*Dt;
end
x=zeros(1,R+1);
x(1,1)=0;
for i=1:R
    x(1,i+1)=x(1,i)+Dt;
end

figure;
plot(x,N,'b',x,P,'r',x,No,'g');