%
% Runge-Kutta integration for a 2D domain 
% Ricardo Carretero Feb 2003
%
function [ru,rv] = rk4_2d_noflux_periodic(u,v,Nx,Ny,dt,oos,dxx,dyy)

 global allt allx ally allu allv;
 global D1 D2 allt allx ally allu allv ran;
 global aa bb gg alpha rho K

%%%%%%%%%%%%%%%%%%%%%%%%%  k1  %%%%%%%%%%%%%%%%%%%%%%%%%
   sdu1 = diff2d_noflux_periodic(u,Nx,Ny,dxx,dyy);
   sdv1 = diff2d_noflux_periodic(v,Nx,Ny,dxx,dyy);
   k1u=dt*(D1*sdu1+reaction_u(u,v));
   k1v=dt*(D2*sdv1+reaction_v(u,v));

%%%%%%%%%%%%%%%%%%%%%%%%%  k2  %%%%%%%%%%%%%%%%%%%%%%%%%
   au=u+0.5*k1u;
   av=v+0.5*k1v;

   sdu1 = diff2d_noflux_periodic(au,Nx,Ny,dxx,dyy);
   sdv1 = diff2d_noflux_periodic(av,Nx,Ny,dxx,dyy);
   k2u=dt*(D1*sdu1+reaction_u(au,av));
   k2v=dt*(D2*sdv1+reaction_v(au,av));

%%%%%%%%%%%%%%%%%%%%%%%%%  k3  %%%%%%%%%%%%%%%%%%%%%%%%%
   au=u+0.5*k2u;
   av=v+0.5*k2v;

   sdu1 = diff2d_noflux_periodic(au,Nx,Ny,dxx,dyy);
   sdv1 = diff2d_noflux_periodic(av,Nx,Ny,dxx,dyy);
   k3u=dt*(D1*sdu1+reaction_u(au,av));
   k3v=dt*(D2*sdv1+reaction_v(au,av));

%%%%%%%%%%%%%%%%%%%%%%%%%  k4  %%%%%%%%%%%%%%%%%%%%%%%%%
   au=u+k3u;
   av=v+k3v;

   sdu1 = diff2d_noflux_periodic(au,Nx,Ny,dxx,dyy);
   sdv1 = diff2d_noflux_periodic(av,Nx,Ny,dxx,dyy);
   k4u=dt*(D1*sdu1+reaction_u(au,av));
   k4v=dt*(D2*sdv1+reaction_v(au,av));

   ru = u + oos*(k1u+2*k2u+2*k3u+k4u);
   rv = v + oos*(k1v+2*k2v+2*k3v+k4v);

