r: right hand side of the Runge-Kutta scheme (-Jacobian+Laplacian terms)
function rhs(nx,ny,dx,dy,re,w,s,r)
aa = 1.0/(re*dx*dx)
bb = 1.0/(re*dy*dy)
gg = 1.0/(4.0*dx*dy)
hh = 1.0/3.0
for i = 2:nx for j = 2:ny