function [x,lamda,s] = mehrotra(A,b,c,x0,lamda0,s0,epsilon,M)

% 1. wähle .... und setze k=0
k=0

% 2. setze (x,lamda,s)=(x_k,lamda_k,s_k), my_k=x'*s/n
x=x0
lamda=lamda0
s=s0
n=size(A,2)
m=size(A,1)
my_k=(x'*s)./n

% 3. falls Norm < epsilon, Norm < epsilon, my< epsilon, Stop

while (~((norm(A*x-b)<epsilon)&&(norm(A'*lamda+s-c)<epsilon)&&(my_k<epsilon)))
    while ((norm(x)<=M)&&(norm(s)<=M))
        my_k=(x'*s)./n
       %5. Löse Gleichung
        % n ist die Größe von A, S,X, I, e
        A
        S=diag(s)
        X=diag(x)
        Z=zeros(n,n)
        I=eye(n,n)
        Z2=zeros(m,m)
        Z3=zeros(m,n)
        Z4=zeros(n,m)
        C=([Z A' I; A Z2 Z3; S Z4 X])
        delta_xN = sym('delta_xN','real')
        delta_lamdaN = sym('delta_lamdaN','real')
        delta_sN = sym('delta_sN','real')
        
        c1=[delta_xN, delta_lamdaN, delta_sN]' % Spaltenvektor
        e=ones(n,1)
        D=[A'*lamda + s - c ; A*x - b ; X*S*e ]
        c1=solve('C(:,(1:n))*delta_xN + C(:,(n+1):(n+m))*delta_lamdaN + C(:,(n+m+1):(n+m+n))*delta_sN = D')
         % Lösung ist leere Menge
        
        %6. Berechne maximal mögliche Schrittweite entlang delta_xN,
        %delta_sN
        for i=1:n
            if delta_xN(1,i)<0
                q1(i)=-(x(1,i)/delta_xN(1,i))
            end
        end
        % führe q als Hilfsvektor ein, als Quotient
        % wähle für alpha_xN nur die Elemente aus q1
        alpha_xN=min(1, min(q1));
        
        % gleiches für alpha_sN:
        
        for i=1:n
            if delta_sN(1,i)<0
                q2(i)=-(s(1,i)/delta_sN(1,i))
            end
        end
        alpha_sN=min(1, min(q2));
        
        % 7. Setze my_n_t= ...
        my_N_t=((x + alpha_xN*delta_xN )'* (s + alpha_sN*delta_sN ))./n
        
        my_c=my_k*((my_N_t/my_k)^3)
        
        % 8. Löse
        % C wie oben, D neu, c neu
        delta_XN=diag(delta_xN);
        delta_SN=diag(delta_sN);
        D=-( [A'*lamda+s-c; A*x-b; X*S*e -my_c*e + delta_XN*delta_SN*e] );
        c2=([delta_x_c, delta_lamda_c, delta_s_c])';
        solve('C*c2=D');
        
        % 9. Wähle Dämpfungsparameter und berechne primale und duale
        % Schittweite entlang delta_xc, delta_sc
        
        ny_k=0.9995 % Darf ich das???
        
        
         for i=1:n
            if delta_x_c(1,i)<0
                q3(i)=-(x(1,i)/delta_xc(1,i))
            end
         end
                    
        alpha_c_max_x= min(q3);
        
         for i=1:n
            if delta_s_c(1,i)<0
                q4(i)=-(s(1,i)/delta_sc(1,i))
            end
         end
        
        alpha_c_max_s= min(q4);
        alpha_c_x=min(1, ny_k*alpha_c_max_x);
        alpha_c_s=min(1, ny_k*alpha_c_max_s);
        
        x  =   x   +  alpha_c_x*delta_x_c;
        lamda=lamda + alpha_c_s*delta_lamda_c;
        s= s +  alpha_c_s*delta_s_c;
        
        k=k+1
      
    end
end