clear all; close all;

L=2*pi; n=64;
x2=linspace(-L/2,L/2,n+1); x=x2(1:n); y=x; 
dx=x(2)-x(1);
[X,Y]=meshgrid(x,y);
u=exp(-X.^2-Y.^2);
v=exp(-X.^2-Y.^2);
u2=reshape(u,n^2,1); v2=reshape(v,n^2,1);
D1=0.1;
uv=[u2; v2];

e1=ones(n,1);
A=spdiags([e1 -2*e1 e1],[-1 0 1],n,n);
A(1,n)=1; A(n,1)=1;

I=eye(n);
B=kron(A,I)+kron(I,A);

[t,uvsol]=ode45('rhs',0:0.1:0.8,uv,[],D1,B,dx,n);

for j=1:length(t)
  subplot(3,3,j)
  pcolor(X,Y,reshape(uvsol(j,1:n^2),n,n));
  shading interp
end