clear;
close;
dt=0.01; %time step
M=1000; %number of discrete space points
N=2000; %total number of time steps
lambda=200; % wavelength of ripple

P=2;
Q=.2;
R=9;

Abar=R+1;
Ibar=(R+1)^2;

Ahat=0.1*Abar;
Ihat=0.1*Ibar;

s=1:M;

A0=Abar*ones(1,M)+Ahat*cos(2*pi*s/lambda);
I0=Ibar*ones(1,M)+Ihat*cos(2*pi*s/lambda);

A=A0;
I=I0;

for t=1:N,
    Aper=[A(M) A A(1)];
    Iper=[I(M) I I(1)];
   
    ddA=diff(diff(Aper));
	ddI=diff(diff(Iper));
	dAdt=1+R*A.^2./I-A+ddA;
	dIdt=Q*(A.^2-I)+P*ddI;

    A=A+dAdt*dt;
    I=I+dIdt*dt;
    subplot(1,2,1)
    plot(s,A,'go-');
    axis([0 1000 0 30]);
    title('Activator');
    xlabel('Postion')
    ylabel('Activator Concentration');
    subplot(1,2,2)
    plot(s,I,'rx-');
    axis([0 1000 0 400]);
    title('Inhibitor');
    xlabel('Postion')
    ylabel('Inhibitor Concentration');
    drawnow;
end;