function [k_ideal,C] = plug_MSE_diff(lambda,b0h,sigma,maxk,stepk)

MSE(1) = inf;

for kk=0:maxk
    
    k = kk*stepk;           
    VAR = (sigma^2)*sum(lambda.^(2*k));    
    BIAS2= sum((b0h.^2).*(1-lambda.^k).^2);   
    MSE(kk+2) = VAR + BIAS2;
    
    if MSE(kk+2) > MSE(kk+1)
        break;
    end
        
end
    
[C,k_diff] = min(MSE(:));

k_ideal = stepk*(k_diff-2);