Mesh flattening/ smoothing by average of points

Question

I am trying to create the above effect of smoothing or flattening a mesh. This mesh boundary is z = 2*exp(-(y-5).^2).*sin(x) + exp(-x.^2).*cos(y). The flattening should be done by creating a new z value of an interior point (excluding boundary points), where the new value is the average of 3x3 grid points centered at the point. I believe the z values of boundary points should not change.

Am I calculating the average incorrectly?

[x,y] = meshgrid(-10:0.5:10,-10:0.5:10);
z = 2*exp(-(y-5).^2).*sin(x) + exp(-x.^2).*cos(y);

while true
    clf
    surf(z);
    ylim([-2,2])
    n = input('Press enter to continue.');
    sx0 = size(x);
    sx = sx0(1);
    sy = size(y);
    ix = sx0(2); % number of elements along the x-axis
    iy = sy(1); % number of elements along the y-axis
    z1 = z;
    % for each interior point
    for i = 1:(ix-1)
        for j = 2:iy
            %compute the average for a 3x3 grid points.
           z1(j, i) = (sum(z((i-1):(i+1),(j-3):(j+3)))); %ERROR HERE
           pause(0.5);
        end
    end
end

Answer 1

In addition to the point about not updating z, your indexing is all messed up. The updating step should be something like this:

    ni = size(z, 1);
    nj = size(z, 2);
    z1 = z;
    % for each interior point
    for i = 2:ni-1
        for j = 2:nj-1
           z1(i, j) = sum(sum(z(i-1:i+1, j-1:j+1))) / 9;
        end
    end
    z = z1;

Answer 2

I've figured out the problem, it was mostly with my z1 declaration. I was trying to keep the previous value of z but ended up totally not updating z. I've added a z0 to fix that.

[x,y] = meshgrid(-10:0.5:10,-10:0.5:10);
z = 2*exp(-(y-5).^2).*sin(x) + exp(-x.^2).*cos(y);
z0 = z;

while true
    clf
    surf(z);
    zlim([-2,2]);
    figure(1);
    n = input('Press enter to continue.');
    sx = size(x)*1;
    sy = size(y);
    ix = sx(2); % get the number of elements along the x-axis
    iy = sy(1); % get the number of elements along the y-axis
    z1 = z;
    for t = 0:30 %0:time in seconds
        for i = 2:(ix-1)
            for j = 2:(iy-1)
                %compute the average for a 3x3 grid points.
                z1(i,j) = (mean(z((i-1):(i+1), (j-1):(j+1)),'all'));
            end
        end
        z = z1;
        surf(z1);
        pause(0.5);
    end
    z=z0;
end