2D Crank-Nicolson ADI scheme

I want to solve the next PDE system using a Crank-Nicolson scheme:

where Du and Dv are diffusive constants, and "a" and "b" are just positive constants. It is possible to find simulations for this PDE system using Crank-Nicolson scheme?. I have the 1D dimensional problem solved and here is the code:

   % Modelo de Sel'kov en 1D
   % Este código en MATLAB simulara el modelo de Sel'kov para la glucólisis en 1D
   clear all
   close all

   %%%%%%%%%%%%%%%%%%%%%
   %%%Inicialización %%%
   %%%%%%%%%%%%%%%%%%%%%

   % Parámetros
   bc  = 0.6; %0.98; %0.90033; %0.9502; % %0.6;
   ac = 0.08; %0.08; %0; %
   Dv = 0.001; %0.02; %0.002; 
   Du = 0.001;  %Constantes difusivas

% Información inicial y mallado
 % w = 10;  %sin patrón
 w = 1;  % patrón

   Nx = 1000; % Total de puntos discretizados en el dominio [0,L]
   x = linspace(0,w,Nx);
   dx = x(2) - x(1);
   %pause();

   dt = 0.01; % tamaño del paso
   t = 0:dt:100;
   Nt = length(t); %numero de puntos en el tiempo

   % Condiciones para la superficie
   [X, T] = meshgrid(x, t);
   %U = 0*X;
   %V = 0*X;

   % Vectores columna (más fáciles)
   x = x(:);
   t = t(:);

   %Condición inicial: pequeña perturbación lejos del estado estable
   u =  bc.*ones(length(x),1) + 0.5.*rand(Nx, 1); %0.01.*(cos(pi.*x)) +    %exp(-(x-0.5).^2);
   v =  (bc/(ac + bc^2)).*ones(length(x),1) + 0.5.*rand(Nx,1); %0.01.*(cos(2.*pi.*x)) +  % exp(-(x-0.5).^2);

   % Condiciones iniciales salvadas.
   U(1,:) = u;
   V(1,:) = v;

   %%%%%%%%%%%%%%%%%%%%%%%%%
   %%%Generando la matriz%%%
   %%%%%%%%%%%%%%%%%%%%%%%%%
  % Usando un esquema implícito 
  %  Para u
    lambda1 = dt.*Du./(2.*dx.^2);
    Au = eye(Nx).*(1+2.*lambda1);
    Au(1+1:Nx+1:end) = -lambda1;
    Au(Nx+1:Nx+1:end) = -lambda1;
    Au([1,end]) = 1+lambda1;
    Bu = eye(Nx).*(1-2.*lambda1);
    Bu(1+1:Nx+1:end) = lambda1;
    Bu(Nx+1:Nx+1:end) = lambda1;
    Bu([1,end]) = 1-lambda1;

%  Para v
  lambda2 = dt.*Dv./(2.*dx.^2);
  Av = eye(Nx).*(1+2.*lambda2);
  Av(1+1:Nx+1:end) = -lambda2;
  Av(Nx+1:Nx+1:end) = -lambda2;
  Av([1,end]) = 1+lambda2;
  Bv = eye(Nx).*(1-2.*lambda2);
  Bv(1+1:Nx+1:end) = lambda2;
  Bv(Nx+1:Nx+1:end) = lambda2;
  Bv([1,end]) = 1-lambda2;

   %%%%%%%%%%%%%%%%%%%%%%%%%

Resultados %%%

   %%%%%%%%%%%%%%%%%%%%%%%%%

   axis([0 1 0 1]) 
   uma(1)=u(Nx/2);
   vma(1)=v(Nx/2);
   for j = 1:Nt-1
       % f y g son los términos no lineales en el modelo de Sel'kov

       f= -u+ac.*v+u.^2.*v;
       g= bc-ac.*v-u.^2.*v;   

       %f = u.^2./v-bc*u;
       %g = u.^2 - v;
       %size(Bu), size(u),
       % En cada paso resolvemos el sistema de ecuaciones
       u = Au\(Bu*u + dt.*f); %u = Bu\(u + dt.*f);   
       v = Av\(Bv*v + dt.*g); %v = Bv\(v + dt.*g);

       uma(j+1)=u(Nx/2);
       vma(j+1)=v(Nx/2);
       %u = up;
       %v = vp;

       % Gráficas
       plot(x,u,'g.-', 'linewidth',1);
       hold on;
       plot(x,v,'r.-', 'linewidth',1);
       hold off;
       legend('ATP','ADP')
       axis([0 1 0 2])
       title(['t = ', num2str(j*dt)],'fontsize',24)
       drawnow;
       %pause();
       % Datos para la superficie
       U(j+1,:) = u;
       V(j+1,:) = v;
   end

   %%%%Gráfica de la superficie %%%%
   figure(1);
   s = surf(x, t, U);
   set(s, 'EdgeColor', 'none', 'FaceColor', 'interp'); 
   xlabel('<----x---->')
   ylabel('<----t---->')
   zlabel('<----u---->')

    figure(2);
   s = surf(x, t, V);
   set(s, 'EdgeColor', 'none', 'FaceColor', 'interp'); 
   xlabel('<----x---->')
   ylabel('<----t---->')
   zlabel('<----v---->')

   %%%%contour plot %%%
   figure(3);
   p = pcolor(x, t, U);
  set(p, 'EdgeColor', 'none', 'FaceColor', 'interp');

figure(4);
q = pcolor(x, t, V);
set(q, 'EdgeColor', 'none', 'FaceColor', 'interp');

figure(5);
 plot(t,uma,'g.-', t,vma,'r.-');
 legend('ATP','ADP')

How I should modify my 1D code in order to solve the 2D problema put at the top of this question?