% spatial domain of x and y 
Lx=20; Ly=20; 

% discretization points in x and y 
nx=100; ny=nx; 

% build a vector of ones  
e1=ones(nx,1); 

% construct spatial discretization matrix A 
A=spdiags([e1 -2*e1 e1],[-1 0 1],nx,ny); 

% periodic boundaries: 
A(1,ny)=1; A(nx,1)=1; 

% idendity matrix of size nx^2 
I=eye(nx);

% 2D differentiation matrix using A from 1D 
B=kron(I,A)+kron(A,I); 

% account for periodicity 
x2=linspace(-Lx/2,Lx/2,nx+1); 
x=x2(1:nx); 

% account for periodicity 
y2=linspace(-Ly/2,Ly/2,ny+1); 
y=y2(1:ny); 
dx = x(2)-x(1); % or: dx = 1/n 

% set up for 2D initial conditions 
[X,Y]=meshgrid(x,y); 

% generate a Gaussian initial condition matrix 
U=exp(-X.^2-Y.^2); 

% elements in reshaped initial condition 
N=nx*ny;

% reshape into a vector 
u=reshape(U,N,1); 

% solve the ode for kappa = 1 
k = 1; 
tspan = [0 1]; 
[t,y]=ode45('rhs2D',tspan,u,[],k,dx,B); 

% plot the solution for the last time
figure(3)
set(gca,'FontSize',20);
fin_sol = y(end,:);
u_final = reshape(fin_sol,nx,nx); 
surf(X,Y,u_final);
xlabel('x'); ylabel('y'); zlabel('u');