Floyd Steinberg Dithering Matlab - What am I doing wrong?

I am trying to implement Floyd Steinberg Dithering in MATLAB, using the pseudocode on the Wikipedia page https://en.wikipedia.org/wiki/Floyd%E2%80%93Steinberg_dithering

My code is below

image = double(imread("dithering.jpg")) ./ 255;
levels = 2;

image_quantised = round(image .* (levels - 1)) ./ (levels - 1);
error = image - image_quantised;

height = size(image(:, :, 1), 1);
width = size(image(:, :, 1), 2);

image_dithered = image_quantised;

for y = 1:height - 1
    for x = 2:width - 1
                
        image_dithered(y    , x + 1, :) = image_dithered(y    , x + 1, :) + error(y, x, :) .* 7 / 16;
        image_dithered(y + 1, x - 1, :) = image_dithered(y + 1, x - 1, :) + error(y, x, :) .* 3 / 16;
        image_dithered(y + 1, x    , :) = image_dithered(y + 1, x    , :) + error(y, x, :) .* 5 / 16;
        image_dithered(y + 1, x + 1, :) = image_dithered(y + 1, x + 1, :) + error(y, x, :) .* 1 / 16;
        
    end
end

imshow(image_dithered)  % Image 1
imshow(dither(mean(image, 3)))  % Image 2

Image 1

Image 2

I am expecting the result in Image 2, but I am getting Image 1. It looks as though the algorithm isn't doing anything. Any ideas? :)

Edit: I have tried initialising image_dithered with different values; all zeros, the quantised image, and the original image. None of them work correctly

Edit 2: I'm getting closer by now calculating the error and quantisation within the loop. Still not spot on however.

for y = 1:height - 1
    for x = 2:width - 1
        new_pixel = round(image_dithered(y, x, :) .* (levels - 1)) ./ (levels - 1);
        image_dithered(y, x, :) = new_pixel;

        error = image(y, x, :) - new_pixel;

        image_dithered(y    , x + 1, :) = image_dithered(y    , x + 1, :) + error .* 7 / 16;
        image_dithered(y + 1, x - 1, :) = image_dithered(y + 1, x - 1, :) + error .* 3 / 16;
        image_dithered(y + 1, x    , :) = image_dithered(y + 1, x    , :) + error .* 5 / 16;
        image_dithered(y + 1, x + 1, :) = image_dithered(y + 1, x + 1, :) + error .* 1 / 16;
    end
end

Edit 3: Thanks for the @saastn and @Cris Luengo, both of the answers helped me work out where I was going wrong and it appears to be working as expected now!

The fixed code is below for completeness.

height = size(image(:, :, 1), 1);
width = size(image(:, :, 1), 2);

image_dithered = image;

for y = 1:height - 1
    for x = 2:width - 1
        old_pixel = image_dithered(y, x, :);
        new_pixel = round(image_dithered(y, x, :) .* (levels - 1)) ./ (levels - 1);
        
        image_dithered(y, x, :) = new_pixel;
        
        error = old_pixel - new_pixel;
        
        image_dithered(y    , x + 1, :) = image_dithered(y    , x + 1, :) + error .* 7 / 16;
        image_dithered(y + 1, x - 1, :) = image_dithered(y + 1, x - 1, :) + error .* 3 / 16;
        image_dithered(y + 1, x    , :) = image_dithered(y + 1, x    , :) + error .* 5 / 16;
        image_dithered(y + 1, x + 1, :) = image_dithered(y + 1, x + 1, :) + error .* 1 / 16;
        
    end
end

imshow(image_dithered)
imshow(dither(mean(image, 3)))

The problem is that you tried to improve the pseudocode and remove oldpixel! Note that the algorithm does not calculate an error between the quantized pixel and its corresponding value in the original image, but rather the error between the quantized pixel and the previous value in the dithered image, which may have been updated while scanning the previous pixels. Bring oldpixel back and review the whole algorithm one more time.

But even after the modifications, you can not expect the results to match the MATLAB output, which could be the result of differences in the details of the two implementations.