function [wel,n] = func_eqm_iter(m)
% Solve the equilibrium by iteration

w_c = 1;        w_u = 1;
r_cc = 1;       r_cu = 1;
r_uu = 1;       r_uc = 1;
X_c = 1;        X_u = 1;
l_cc = 0.5;     l_cu = 0.5;
l_uu = 0.5;     l_uc = 0.5;

diff = 1;
cc = 0;

while diff > m.tol && cc < 100000
    [wel,n,w_c_out,w_u_out,r_cc_out,r_cu_out,r_uu_out,r_uc_out,X_c_out,X_u_out,l_cc_out,l_uu_out,l_cu_out,l_uc_out] ...
        = func_eqm_update(m,w_c,w_u,r_cc,r_cu,r_uu,r_uc,X_c,X_u,l_cc,l_uu,l_cu,l_uc);

    eqm_in  = [w_c;w_u;r_cc;r_cu;r_uu;r_uc;X_c;X_u;l_cc;l_uu;l_cu;l_uc];
    eqm_out = [w_c_out;w_u_out;r_cc_out;r_cu_out;r_uu_out;r_uc_out;X_c_out;X_u_out;l_cc_out;l_uu_out;l_cu_out;l_uc_out];
    diff = max(abs(eqm_out - eqm_in));
    cc   = cc+1;
    % disp(diff)
    % disp(cc)
    
    w_c = m.weight*w_c_out + (1-m.weight)*w_c;
    w_u = m.weight*w_u_out + (1-m.weight)*w_u;
    r_cc = m.weight*r_cc_out + (1-m.weight)*r_cc;
    r_cu = m.weight*r_cu_out + (1-m.weight)*r_cu;
    r_uu = m.weight*r_uu_out + (1-m.weight)*r_uu;
    r_uc = m.weight*r_uc_out + (1-m.weight)*r_uc;
    X_c = m.weight*X_c_out + (1-m.weight)*X_c;
    X_u = m.weight*X_u_out + (1-m.weight)*X_u;
    l_cc = m.weight*l_cc_out + (1-m.weight)*l_cc;
    l_uu = m.weight*l_uu_out + (1-m.weight)*l_uu;
    l_cu = m.weight*l_cu_out + (1-m.weight)*l_cu;
    l_uc = m.weight*l_uc_out + (1-m.weight)*l_uc;

end

end

