%                     Simulation of Reaction Diffusion System
clear
% Parameter
Diff=1e-4;
% Numerical Parameters
n=101;
dt=0.01;
tmax=6;
% Plotting Parameters
tplot=0.1;

% SET UP initial values
x=linspace(0,1,n);
u=0.5*sin(x)';

% SET UP Figures
clf   %  clear figures
pict=plot(x,u,'b','EraseMode','xor');
axis([0 1 0 1.2]);

% step in time
t=0;
nsteps=0;
tp=tplot;
while t<tmax
   [u]=stepRD(u,dt,Diff);
   t=t+dt;
  %storage
   nsteps=nsteps+1;
   bigu(:,nsteps)=u; % could speed up by allocating matrix before loop
   bigt(nsteps)=t;
  %plots
   if t>=tp 
       tp=tp+tplot;
       set(pict,'Xdata',x,'Ydata',u)
       drawnow;
   end
end
disp('Simulation ended')
%  Try plotting u vs t for the point x=0.87 (or so)
% at the command line, enter:
%   clf
%   plot(bigt,bigu(87,:))
%  Try plotting u vs x each unit of time (t=0,1,2,3,4,5,6)
% at the command line, enter:
%   clf
%   plot(x,bigu(:,1:100:end))  % works for dt=0.01 so 100 timesteps => t=1
