function div = divergence(u,v,dx,dy)
[n,p] = size(u);
dudx=u;dudx(:,:)=0;
dvdy=dudx;

u(:,end)=0;
dudx(:,2:end)=(u(:,2:end)-u(:,1:end-1))/dx;
dudx(:,1)=u(:,1)/dx;

v(end,:)=0;
dvdy(2:end,:)=(v(2:end,:)-v(1:end-1,:))/dy;
dvdy(1,:)=v(1,:)/dy;

div=dudx+dvdy;
