V=1*10^10;
Q=16438356.16;
Dt=1;  %Dt
R=300
T=V/Q;     %XRONOS PARAMONIS
L=1.5*10^8; %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)=21.42;    %ug/L
P(1,1)=2.73;     %ug/L
No(1,1)=2.38;    %ug/L
Kg=0.04;         %L/ug/day
Ko=0.48;      %1/day
K1=0.60;       %1/day
Ks=0.26;       %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');