Multithresholding using Otsu method with 3 clusters, not the same threshold with own code

Hello,

I wanted to use the Matlab function "multithresh" to get a correct image segmentation. Next to it, I created my own thresholding code using the same (Otsu) method. The problem is that the thresholds obtained are not the same. I could not find an error in my code. My code is written rather mathematically and is time consuming but seems, compared to a global segmentation, to yield similar results.

I'm now in doubt about the results of my code and the multithresh. Would it be possible to have a look at my code and compare the results?

Thank you for your precious help, Olivier

 % Amount of pixels for each intensity
 for i = 1:size(image,1)
    for j = 1:size(image,2)
        for k = 1:size(L)
            if image(i,j) == L(k)
                n(k)=n(k)+1;
            end
        end
    end
 end

 %Total amount of pixels for all intensities
 MN = size(image,1)*size(image,2);

 %Normalized amount of pixels for each intensity
 p = n/MN;

 %Compute mean global intensity
 mG = 0;
 for i = 1:256
    mG = mG + L(i)*p(i);
 end

 %Global variance
 varG = 0;
 for i = 1:256
    varG = varG + (L(i)-mG)^2*p(i);
 end

 %Loop to find maximal in-between variance
 %Loop with condition: 0 < k1 < k2 < 255 (corresponds to nuances and not positions)
 varB = zeros(256,256);
 for a = 3:255              % position
    k2 = a-1;               % nuances
    for b = 2:k2            % position
        k1 = b-1;           % nuances

        %Probability of classes Pr1, Pr2, Pr3
        Pr1 = 0;
        for i = 1:(k1+1)
           Pr1 = Pr1+p(i);
        end
        Pr2 = 0;
        for i = (k1+2):(k2+1)
           Pr2 = Pr2+p(i);
        end
        Pr3 = 0;
        for i = (k2+2):256
           Pr3 = Pr3+p(i);
        end    

        %Compute mean intensity of each class
        m1 = 0;
        for i = 1:(k1+1)
           m1 = m1 + 1/Pr1*L(i)*p(i);
        end
        m2 = 0;
        for i = (k1+2):(k2+1)
           m2 = m2 + 1/Pr2*L(i)*p(i);
        end
        m3 = 0;
        for i = (k2+2):256
           m3 = m3 + 1/Pr3*L(i)*p(i);
        end

        %In between-class variance of actual k1 and k2
        varB(k1,k2) = Pr1*(m1-mG)^2+Pr2*(m2-mG)^2+Pr3*(m3-mG)^2;

    end
 end

 %Maximal in between-class variance 
 varB_max = max(max(varB));