Исправить ошибку в коде по приближению значений

%дано: 
t=0.1; 
a=10; 
d=0.1;
% Начальные приближения искомых параметров системы
    k0=1;
    k1=1;
    Mx=1;
    Dx=1;
    Mxmin=1;
    ind=1;
    nomer=0;
    for k=4:0.1:20
        p=1;
        while p>0.001
 
             k01=(2/mx)*erf(real((mx/((dx)^0.5))));
            k11=(1/(2*(dx^0.5)))*(sqrt(1-4*(erf(real(mx/((dx)^0.5))))^2)+(2/sqrt(2*pi))*exp(-(1/2)*(mx/(dx^0.5))^2));
            mx1=2/(k0*k);
            dx1=((k*k1*d*(t*a+1))/(t*(a^2)+a+k1*k))^0.5;
 
            pk0= (abs(k01-k0))^2;
            pk1=(abs(k11-k1))^2;
            pMx= (abs(Mx1-Mx))^2;
            pDx=(abs(Dx1 - Dx))^2;
            p2=pk0+pk1+pMx+pDx;
            if p2<p
               p=p2;
            end;   
            Mx=real(Mx1);
            Dx=real(Dx1);
            k0=k01;
            k1=k11;
        end;    
        
      % Запись результатов в массивы
        kmas(ind)=k;
        Mxmas(ind)=Mx;
        Dxmas(ind)=Dx;
        s1=num2str(Mx);
        s2=num2str(sqrt(Dx));
        s =strcat('normcdf(M,',s1,',',s2,')-normcdf(-M,',s1,',',s2,')-0.8');
        M(ind)=fsolve(inline(s),0);
       if Mxmin>M(ind)
          Mxmin=M(ind);
          nomer=ind;
       end;
        ind=ind+1;
  end;      
 
 %построение графика зависимости мат.ожидания от k
 figure
plot(kmas,Mxmas)
grid on;
Xlabel('k');
Ylabel('Mx');
 
 %построение графика зависимости дисперсии от k
 figure
plot(kmas,Dxmas)
grid on;
Xlabel('k');
Ylabel('Dx');
 
 %построение графика зависимости M от k
 figure
plot(kmas,M)
grid on;
Xlabel('k');
Ylabel('M');
 
%вывод данных
Mxmin
Dx  
k=kmas(nomer)

Undefined function or variable 'mx'.

% ИСХОДНЫЕ ДАННЫЕ СОГЛАСНО ВАРИАНТУ
a=5; 
d=1; 
T1=0.5;
T2=3;
% ЗАДАЁМ С 0 ДО 1
    k0=1;
    k1=1;
    Mx=1;
    Dx=1;
    M_min=1;
    i=1;
    nom=0;
    for k=4:0.001:20% НАЧИНАЕМ С 4, ТАК КАК ТОЧНОСТЬ МЕТОДА НЕ ПОЗВОЛЯЕТ ВЫБРАТЬ БОЛЬШИЙ ДИАПАЗОН С ЗАДАННЫМ ШАГОМ   
        N=1;% ЗАДАЁМ ЗНАЧЕНИЕ НЕВЯЗКИ
        while N>0.01
          
            % СИСТЕМА НЕЛИНЕЙНЫХ УРАВНЕНИЙ (29)
            k01=(2/Mx)*erf(real((Mx/((Dx)^0.5))));
            k11=(1/(2*(Dx^0.5)))*(sqrt(1-4*(erf(real(Mx/((Dx)^0.5))))^2)+(2/sqrt(2*pi))*exp(-(1/2)*(Mx/(Dx^0.5))^2));
            Mx1=2.5/(1+k0*k);
            Dx1=((k1*k1*k*k*d*(T1*T2*(T1+T2+T1*T2*a)))/(T1*T2*(1+k*k1)*((T1+T2+T1*T2*a)*(T1*a+T2*a+1+k*k1)+(T1*T2*a+k*k1*a*T1*T2))))^0.5;
            
            N1=(k01-k0)^2+(k11-k1)^2+(Mx1-Mx)^2+(Dx1-Dx)^2;
            if N1<N
                N=N1;
            end;
           Mx=real(Mx1);
           Dx=real(Dx1);
           k0=k01;
           k1=k11;
        end;
        
    % ФОРМИРОВАНИЕ МАССИВОВ
    % Преобразование массива в строку символов st1-2 с точностью до четырех десятичных разрядов и экспоненциальным представлением
    % Горизонтальное объединение соответствующих рядов массивов символов
    % st1, st2 в единный массив str(ряд) результирующего массива. В normcdf передаём неизвестное М
    % st1 - массив(строка) значений мат.ожидания от k
    % st2 - массив(строка) значений СКО от k
        kmas(i)=k;
        M_mas(i)=Mx;
        Dxmas(i)=Dx;
        st1=num2str(Mx);
        st2=num2str(sqrt(Dx));
        st=strcat('normcdf(M,',st1,',',st2,')-normcdf(-M,',st1,',',st2,')-0.8');
        M(i)=fsolve(inline(st),0);
       if M_min>M(i)
           M_min=M(i);
           nom=i;
       end;
        i= i+1;
end;
 
% Построение графика зависимости М от k
figure
plot(kmas, M, '.c')
grid on;
title('M(k)'); 
xlabel('k'); 
ylabel('M'); 
% Построение графика зависимости дисперсии от k
figure
plot(kmas,Dxmas, '.r')
grid on;
title('Dx(k)'); 
xlabel('k'); 
ylabel('Dx'); 
% Построение графика зависимости Мат.ожидание от k
figure
plot(kmas, M_mas, '.b')
grid on;
title('Mx(k)'); 
xlabel('k'); 
ylabel('Mx'); 
% Вывод данных:
Mx=real(Mx)
Dx=real(Dx)
M=real(M_min)
k=kmas(nom)