im studying travelling waves for multidimension reaction diffusion equation so i start with 2D Fisher KPP to do , i use implicit Crank-Nicolson Finite diffrence sheme

here the code:

r = 1.0;
D = 1.0;

% x in [x0,x1]
h = 0.025; 
Sx = 10.0;
Sy = 6.0;
x0 = 0;  
x1 = Sx;
y0 = 0;
y1 = Sy;
x = x0:h:x1; % discretisation
y = y0:h:y1; 
[X,Y] = meshgrid(x,y);
J1 = length(x); 
J2 = length(y);
J  = J1*J2; 

u = zeros(J,1); 
newu = zeros(J,1);

% discretisationof laplacian with Newmann condition on boundary  
% no flux, i.e. du/dx = 0 en x0 et en x1



cornertopleft = 1;
cornerbottomleft = J2;
cornertopright = (J1-1)*J2+1;
cornerbottomright = J;
sideleft = 2:J2-1;
sidetop = J2+1:J2:J2*(J1-2)+1;
sidebottom = 2*J2:J2:J2*(J1-1);
sideright = J2*(J1-1)+2:J-1;
side = [cornertopleft, cornertopright, cornerbottomleft, cornerbottomright, ...
    sideleft, sidetop, sidebottom, sideright];
interieur = setdiff(1:J, side);

% interieur
L0 = sparse(interieur,interieur,-4,J,J); 
L0 = L0 + sparse(interieur,interieur+1,1,J,J);
L0 = L0 + sparse(interieur,interieur-1,1,J,J);
L0 = L0 + sparse(interieur,interieur+J2,1,J,J);
L0 = L0 + sparse(interieur,interieur-J2,1,J,J);

% initial condition 
% u( X(:).^2 + Y(:).^2 < 3 ) = 0.8;  
u = exp( - X(:).^2/0.5 - Y(:).^2/2 ); %smooth initial ondition


% time , simulation parameter 
t0 = 0;
tfinal = 35; 
t = t0;
tp = t0;
k=0.5;
figure(1); clf;
image(x,y,255*reshape(u,J2,J1));
disp('appuyer sur une touche pour continuer');
pause

% Schema implicite Crank-Nicolson
A = (speye(J) - dt/h^2*D/2*L0);

% PRINCIPALE loop
tic
while t < tfinal
    drawnow;
    newu =A(u + k*r*u.*(1-u) + k/h^2*D/2*L0*u);
    newu(sideleft) = newu(sideleft+J2);
    newu(sideright) = newu(sideright-J2);
    newu(sidetop) = newu(sidetop+1);
    newu(sidebottom) = newu(sidebottom-1);
    u = newu;
    if ( fix(10*t) > fix(10*tp) )
        image(x,y,255*reshape(u,J2,J1));
    end
    tp = t;
    t = t + k;
    fprintf("t = %.5fn",t);
end

The code run correcltly however i want to calculate the wave speed ,any suggestions ? i tried to integrate the solution and divide by final time but it didnt give the expeted value so i think that this way is wrong Im expecting a value c >=2