%
% Ricardo Carretero 20/Oct/2003
% Integrates a SYSTEM of 2 reaction diffusion eqs in 1D using finite differences
%
%        u_t = f(u,v) + D1 u_xx
%        v_t = g(u,v) + D2 v_xx
%
% Ex:  PDE_2pop(N,dt,tf,L,R,disp,savesnaps)
%
% npts=N, dt:time step, tf:t_final, domain:[L,R]
%
%Ex: 
%global aa bb dx dxx D1 D2 oos allt allx allu allv;
%N=75;dt=0.002;tf=200;L=0;R=100;snaps=100;disp=100;
%aa=1;bb=3.65;D1=1;D2=335;
%aa=1;bb=3.65;D1=1;D2=345;
%PDE_2pop(N,dt,tf,L,R,disp,snaps)
%wave_plotsurf(allt,allx,allu',R,1,N+1,'no')
%wave_plotsurf(allt,allx,allv',R,1,N+1,'no')
%

function [] = PDE_2pop(N,dt,tf,L,R,exposh,saveframes)

 global allt allx allu allv;
 global dx dxx D1 D2 oos u0 v0;
 global aa bb;
 
 ran=1;
 noise=0;

 maxstep=tf/dt;
 dispscreen= fix(maxstep/10); 

 if(exposh>0) stopdisp= fix(maxstep/exposh); 
 else stopdisp=9e+99;
 end
 if(saveframes>0)stopsave= fix(maxstep/saveframes);
 else stopsave=9e+99;
 end
 
 % Spatial computational domain defined

 dx=(R-L)/N;
 x=L+(0:N-1)*dx;
 dxx = 1/dx^2;
 oos=1/6;
 
 rand('state',0)       %this fixes the seed for the random generation
 
 strlength=0;

 allt = 0:stopsave*dt:tf;
 allu = zeros(saveframes+1,N);
 allv = zeros(saveframes+1,N);
 allx = x;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% display parameters
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 fprintf('\nL \t\t= %f\n',L);
 fprintf('R \t\t= %f\n',R);
 fprintf('dt \t\t= %f\n',dt);
 fprintf('dx \t\t= %f\n',dx);
 fprintf('N \t\t= %f\n',N);
 fprintf('D1 \t\t= %f\n',D1);
 fprintf('D2 \t\t= %f\n',D2);
 fprintf('aa \t\t= %f\n',aa);
 fprintf('bb \t\t= %f\n',bb);
 fprintf('tf \t\t= %f\n',tf);
 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% GENERATING ICs
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 u0 = aa+bb;
 v0 = bb/(aa+bb)^2;
 mode=16;
 bell = exp(-0.00001*(x-(R-L-2)/2).^2);
 uini=u0+  0.1*u0*bell.*(0.0*rand(1,N)+cos(mode*pi*(x-L)/(R-L-dx)));
 vini=v0-0.001*u0*bell.*(0.0*rand(1,N)+cos(mode*pi*(x-L)/(R-L-dx)));
 uini=u0+0.1*u0*(rand(1,N)-0.45);
 vini=v0+0.01*v0*(rand(1,N)-0.5);
 save=0;
 t=0;
 u=uini;
 v=vini;

% %No flux BCs
% u(1)=u(2);
% u(N)=u(N-1);
% v(1)=v(2);
% v(N)=v(N-1);

 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Plot IC
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 figure(1); clf;
 subplot(1,2,1)
 set(gca,'FontSize',[16]);
 plot(x,u,'LineWidth',3); 
 title('u(x,0)');
 subplot(1,2,2)
 set(gca,'FontSize',[16]);
 plot(x,v,'LineWidth',3); 
 title('v(x,0)');
 pause
 
 save=save+1;
 allu(save,:) = u;
 allv(save,:) = v;
     
 fprintf('\nt=%5.2f ',t);
 
 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Main loop
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 for kk=1:maxstep
    
  t=t+dt;
  uN = u+dt*deriv_u(u,v,N);
  vN = v+dt*deriv_v(u,v,N);

%  %Fixed BCs
%  uN(1)=u0;
%  uN(N)=u0;
%  vN(1)=v0;
%  vN(N)=v0;

%  %No flux BCs
%  uN(1)=uN(2);
%  uN(N)=uN(N-1);
%  vN(1)=vN(2);
%  vN(N)=vN(N-1);

  u = uN;
  v = vN;
  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% include noise
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   if(noise>0)
    u  = u+noise*(rand(N)-0.5);
   end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% display time on screen
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   if (fix(kk/dispscreen)==kk/dispscreen)
    fprintf('%5.2f ',t);
    strlength=strlength+1;
    if strlength>=10
     strlength=0;
     fprintf('\n');
    end
   end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% store solution
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   if (fix(kk/stopsave)==kk/stopsave)
    save=save+1;
    allu(save,:) = u;
    allv(save,:) = v;
   end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Plot data
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   if(fix(kk/stopdisp)==kk/stopdisp&kk>1)
     figure(1); clf;
     mytext=['T=(' num2str(kk*dt) ')/(' num2str(tf) ')'];
     subplot(1,2,1)
     set(gca,'FontSize',[16]);
     plot(x,u,'LineWidth',3); 
     text(R-0.2*(R-L),min(u)+(max(u)-min(u))*1.1,mytext);
     subplot(1,2,2)
     set(gca,'FontSize',[16]);
     plot(x,v,'LineWidth',3); 
     text(R-0.2*(R-L),min(v)+(max(v)-min(v))*1.1,mytext);
     pause(0.0001)
   end

 end

 fprintf('\n');
 return;
 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%function rkZn = ode_rk(rkZ)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% global dx dxx oos u0 v0;
% global aa bb;
%
% k1 = dt*deriv_u(rkZ       );
% k2 = dt*deriv_u(rkZ+0.5*k1);
% k3 = dt*deriv_u(rkZ+0.5*k2);
% k4 = dt*deriv_u(rkZ+    k3);
% 
% rkZn = rkZ + oos*(k1 + 2*k2 + 2*k3 + k4);
% 
% return;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function du = deriv_u(u,v,N)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 global dx dxx D1 D2 oos u0 v0;
 global aa bb;
 
%% Periodic BCs:
% up = [ u(1) u(1:N-1) ];
% um = [ u(2:N) u(N) ];

%% Fixed BCs
% up = [ u0 u(1:N-1) ];
% um = [ u(2:N) u0 ];

% No flux BCs
 up = [ u(1) u(1:N-1) ];
 um = [ u(2:N) u(N) ];

 du = f(u,v) + D1*dxx*(up-2*u+um);
 
 return;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function dv = deriv_v(u,v,N)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 global dx dxx D1 D2 oos u0 v0;
 global aa bb;
 
%% Periodic BCs:
% vp = [ v(1) v(1:N-1) ];
% vm = [ v(2:N) v(N) ];

%% Fixed BCs
% vp = [ v0 v(1:N-1) ];
% vm = [ v(2:N) v0 ];

% No flux BCs
 vp = [ v(1) v(1:N-1) ];
 vm = [ v(2:N) v(N) ];

 dv = g(u,v) + D2*dxx*(vp-2*v+vm);
 
 return;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function fuv = f(u,v)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 global dx dxx D1 D2 oos u0 v0;
 global aa bb;

 fuv = aa-u+u.^2.*v;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function guv = g(u,v)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 global dx dxx D1 D2 oos u0 v0;
 global aa bb;

 guv = bb-u.^2.*v;

