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


%%%%%%%%%%%%%%%%%%%%%%%%%  k1  %%%%%%%%%%%%%%%%%%%%%%%%%
   sdx1 = diff2d_fix(u,Nx,Ny,dxx,dyy);
   k1x=dt*(sdx1+reaction(u));

%%%%%%%%%%%%%%%%%%%%%%%%%  k2  %%%%%%%%%%%%%%%%%%%%%%%%%
   a=u+0.5*k1x;

   sdx1 = diff2d_fix(a,Nx,Ny,dxx,dyy);
   k2x=dt*(sdx1+reaction(a));

%%%%%%%%%%%%%%%%%%%%%%%%%  k3  %%%%%%%%%%%%%%%%%%%%%%%%%
   a=u+0.5*k2x;

   sdx1 = diff2d_fix(a,Nx,Ny,dxx,dyy);
   k3x=dt*(sdx1+reaction(a));

%%%%%%%%%%%%%%%%%%%%%%%%%  k4  %%%%%%%%%%%%%%%%%%%%%%%%%
   a=u+k3x;

   sdx1 = diff2d_fix(a,Nx,Ny,dxx,dyy);
   k4x=dt*(sdx1+reaction(a));

   ru = u + oos*(k1x+2*k2x+2*k3x+k4x);

