Solve 2D diffusion equation - ADI Method.

I am trying to solve a 2 dimensional diffusion equation of the form d^2C/dt^t = a(d^2C/dt^x + d^2C/dt^y).

This is to solve an electrochemical problem, where some reactant is cconsumed at one of my boundarues (the electrode). There the concentration is always 0 (Dirichlet Boundary Condition). The opposity boundry in the same direction is a wall (no flux, Neumann Boundary condition). The other two boundaries are at a fixed concentration (Drichlet).

I am discretising this using a finite difference method based on the central difference for the second derivatives, and a backward difference for the first time derivative.

I have the following problems:

1) althoigh I impose symmetrical boundary condistions, the result is highly non-simmerical.

2) The method is supposed to be unconditionally stable, but the stability is limited to values of kxw and kyw < 0.5 , whihc is a feature of exclicit numerical mehods.

My code is below. Any help would be immensely appreciated.

clc; clear all;

% Some initial parameters 
C_O2 =  0.27; %mM
Codml = C_O2/C_O2; %Dimensionless concentration
C_PDMS = 1.7; %mM
Cpdl = C_PDMS/C_O2; %Dimensionless concentration  

D_O2 = 2.4e-9; %m^2s^-1 Difffusuion coefficient 
D_PDMS = 2.1e-9; %m2s-1

Dsimo = C_O2/C_O2;
Dsimp =  D_PDMS/D_O2;

% Dimensions of my domain
Ly =10;
Lx =10; 
T =500;
%Step size chose such that kyw and kxw satisfy the stability criterion (i.e. <0.5)
Dt= 0.2; 
Dx= 0.5;
Dy = 0.5;

%Calculate number of steps i time and mesh points in the X and Y directions 
Nt =round(T/Dt);
Nx = round(Lx/Dx);
Ny = round(Ly/Dy);

% Condition for stabilityof the method: both kyw and kxw need to be smaller than 0.5
kyw = Dsimo*Dt/(2*Dy*Dy);
kxw = Dsimo*Dt/(2*Dx*Dx);

%Pre defining the matrix of solution for the full steps and half steps
MConc = zeros(Ny,Nx,Nt);
CHalf = zeros(Ny,Nx,Nt);

%Initial BC's 
b0  = 1;
MConc(:,:,1) =  b0 + zeros(Ny,Nx);
CHalf(:,:,1) =  b0 + zeros(Ny,Nx);


% Codml =2;
%Dirichlet BC's is Y  (no axial smmety) 
bxl = 0; %left boundary at x =0
bxr = 0; % right boundary at x =  Nx
by0 = 0; %Electrode surafce i.e. bottom boundary at y =0 
byL = 1; % top surface boundary at y = Ny

%fill the matrix with these boundary conditions 
MConc(end,:,:) = byL + zeros(Ny,1,Nt);
MConc(1,:,:) = by0 + zeros(Ny,1,Nt);
MConc(:,end,:) = bxl + zeros(1,Nx,Nt);
MConc(:,1,:) = bxr + zeros(1,Nx,Nt); 

CHalf(end,:,:) = byL + zeros(Ny,1,Nt);
CHalf(1,:,:) = by0 + zeros(Ny,1,Nt);
CHalf(:,end,:) = bxl + zeros(1,Nx,Nt);
CHalf(:,1,:) = bxr + zeros(1,Nx,Nt);

% Creating two tridiagonal matrices, to be used on the half time step and full time step respectively
vecy = ones(Ny-2,1);
vecx = ones(Nx-2,1);

diagMLY = [-kyw (1+2*kyw) -kyw].*vecy;
diagMLX = [-kxw (1+2*kxw) -kxw].*vecx; 

MLY = spdiags(diagMLY,-1:1,Ny-2,Ny-2);
MLX = spdiags(diagMLX,-1:1,Nx-2,Nx-2); 

%Vecotrs containing the Boundary conditions. 
Ry = zeros(Ny-2,1);
Ry(1) = -kyw*by0;
Ry(end) = -kyw*byL;

Rx = zeros(Nx-2,1);
Rx(1) = kxw*bxl;
Rx(end) =  kxw*bxr;

%Loop througn time
for t=2:Nt
     
for yf = 2:Ny-1  
    
   for xf = 2 : Nx-1 
    
       Vy(xf-1,1) =  kyw*MConc(xf,yf-1,t-1) +(1-2*kyw)*MConc(xf,yf,t-1) + kyw*MConc(xf,yf+1,t-1);
    
   end
    
     CHalf(xh,2:Ny-1,t) =  MLY\(Vx - Ry); % Still inside the x loop  Sign or R is minuss because the signs of its elemtns are negative
end 
for xh = 2: Nx-1  
    for yh = 2: Ny-1 
   
    Vx(yh-1,1) =  kxw*MConc(xh-1,yh,t-1) +(1-2*kxw)*MConc(xh,yh,t-1) + kxw*MConc(xh+1,yh,t-1);

    end 

    CHalf(xh,2:Ny-1,t) =  MLY\(Vx - Ry); % Still inside the x loop  Sign or R is minuss because the signs of its elemtns are negative

    
for yf = 2:Ny-1  
    
   for xf = 2 : Nx-1 
    
       Vy(xf-1,1) =  kyw*CHalf(xf,yf-1,t-1) +(1-2*kyw)*CHalf(xf,yf,t-1) + kyw*CHalf(xf,yf+1,t-1);
    
   end
      
   MConc(2:Nx-1,yf,t) = MLX\(Vy + Rx);

end
    
 end  
  
end