function [u,v,phi] = UseCase_Steam_v2(u0,v0,phi0)

    %% Parameter
    M       = 400;
    N       = 250;
    K       = 3000;
    K0      = 3000;
    dK      = 20;  % nach jedem m-ten Zeitschritt wird ein Bild ausgegeben
    tol     = 0.5e-1;    
    direkt  = true;
    Daten   = true;
    
    
    
    fpat = 'Lauf7';
    
    %% Simulationsparameter
    ny      = 0.0;
    dx      = 1;
    dt      = 0.37;
    u_in    = 1;
    v_in    = 0;
    
    alpha = 0.000001;
    
    
    %% Einbau
    Einbau = false(M,N);

    % kleiner zentraler Einbau (Für Strahlenspaltung)
%     x = 75;
%     y = 125;
%     r = 10;
%     for i=x-r:x+r
%         for j=y-r:y+r
%             Einbau(i,j) = ((i-x)^2 + (j-y)^2 <= r^2);
%         end
%     end

    % Für Karman1 mit dt = 0.3
%     x = 60;
%     y = 65;
%     r = 40;
%     for i=x-r:x+r
%         for j=y-r:y+r
%             Einbau(i,j) = ((i-x)^2 + (j-y)^2 <= r^2);
%         end
%     end
% 
%     x = 100;
%     y = 175;
%     r = 25;
%     for i=x-r:x+r
%         for j=y-r:y+r
%             Einbau(i,j) = ((i-x)^2 + (j-y)^2 <= r^2);
%         end
%     end
%     
%     x = 175;
%     y = 125;
%     r = 40;
%     for i=x-r:x+r
%         for j=y-r:y+r
%             Einbau(i,j) = ((i-x)^2 + (j-y)^2 <= r^2);
%         end
%     end




    % 7 kleine Einbauten
%     x = 150;
%     y = 125;
%     r = 20;
%     for i=x-r:x+r
%         for j=y-r:y+r
%             Einbau(i,j) = ((i-x)^2 + (j-y)^2 <= r^2);
%         end
%     end
%     
%     x = 250;
%     y = 125;
%     r = 20;
%     for i=x-r:x+r
%         for j=y-r:y+r
%             Einbau(i,j) = ((i-x)^2 + (j-y)^2 <= r^2);
%         end
%     end
%     
%     x = 75;
%     y = 70;
%     r = 20;
%     for i=x-r:x+r
%         for j=y-r:y+r
%             Einbau(i,j) = ((i-x)^2 + (j-y)^2 <= r^2);
%         end
%     end
%     
%     x = 175;
%     y = 70;
%     r = 20;
%     for i=x-r:x+r
%         for j=y-r:y+r
%             Einbau(i,j) = ((i-x)^2 + (j-y)^2 <= r^2);
%         end
%     end
%     
%     x = 75;
%     y = 175;
%     r = 20;
%     for i=x-r:x+r
%         for j=y-r:y+r
%             Einbau(i,j) = ((i-x)^2 + (j-y)^2 <= r^2);
%         end
%     end
%     
%     x = 175;
%     y = 175;
%     r = 20;
%     for i=x-r:x+r
%         for j=y-r:y+r
%             Einbau(i,j) = ((i-x)^2 + (j-y)^2 <= r^2);
%         end
%     end

    % Trichter
    x = 110;
    y = 115;
    b = 25;
    for i=1:40
        Einbau(x-i,y-i) = true;
        Einbau(x-i,y+b+i) = true;
    end
    
    %% Visualisierung: Rauch-Drüsen
    Duese1 = false(M,N);
    Duese2 = false(M,N);
    b = 10;
    d = 40;
     pt_load = 20;
     pt_empt = 100;
%     Duese(1,100:120) = true(1,21);
    % links sind Düsen
    for i=1:250
        if(mod(i,d)==0 )
            Duese2(1,i:i+b-1) = true(1,b);
        end
    end
    Duese2 = Duese2(1:M,1:N);
    Duese1(1,:) = ones(1,N);
    
    %% Hilfsmatrizen
    Null = zeros(M,N);
    Eins = ones(M,N);
    
    
    %% Anfangsverteilung
    if(nargin == 3)
        u = u0;
        v = v0;
        phi = phi0;
    else
        %% Init
         u      = 0*      ones(M,N);
         v      = 0*      ones(M,N);   
         phi    = 0*      ones(M,N);
    end
    
    
    %% Strömungscharakteristika
    ny_til  = dx^2/dt + ny;
    L       = min(M,N) * dx; 
    
    Re = sqrt(v_in^2 + u_in^2) * L / ny_til
    
    
    %% Randbedingungen
    % Matrix mit 1-Einträgen für Ränder
    S           = zeros(M,N);
    S(1,:)      = ones(1,N);
    S(M,:)      = ones(1,N);
    S(:,1)      = ones(M,1);
    S(:,N)      = ones(M,1);
    S(Einbau)   = Eins(Einbau);
    % Einstrom-Flächen
    u(1,:)      = 0 * ones(1,N);
    u(M,:)      = 0 * ones(1,N);
    v(1,:)      = 0 * ones(1,N);
    v(M,:)      = 0 * ones(1,N);
    % No-Slip-Ränder
    u(:,1)      = 0 * ones(M,1);
    u(:,N)      = 0 * ones(M,1);
    v(:,1)      = 0 * ones(M,1);
    v(:,N)      = 0 * ones(M,1);
    % Einbau
    u(Einbau)      = Null(Einbau);
    v(Einbau)      = Null(Einbau);
    
    % Die Randboolwerte  werden in einen Vektor überführt
    for i=1:M
        vS((i-1)*N+1:i*N,1)         = S(i,:)';
    end
    for i=1:M
        vS(M*N+(i-1)*N+1:M*N+i*N,1) = S(i,:)';
    end
    vNS = ones(2*M*N,1) - vS;
    vNS = (vNS == 1);
    
  
    if(false)
        %% Anfangsverteilung plotten
        figure;
        clf;
        hold on;
        grid off;
        axis equal;
        title('Anfangsverteilung');
        xlabel('x');
        ylabel('y');
        colormap((jet+white)/2);
        pcolor(((u .* u + v .* v).^0.5)');
        shading interp;
        quiver(u',v','black');
    end
    
    if(false)
        %% Für Detail-Ansichten
        m1 = 1
        m2 = 400
        m  = 100
        n1 = 1
        n2 = 250
        n  = 62

        figure;

        colormap((jet+white)/2);
        hold on;
        grid off;
        axis equal;
        title('Anfangsverteilung');
        xlabel('x');
        ylabel('y');

        [X,Y] = meshgrid(m1:ceil((m2-m1)/m):m2, n1:ceil((n2-n1)/n):n2);
        norm_u = ((u(m1:ceil((m2-m1)/m):m2,n1:ceil((n2-n1)/n):n2) .* u(m1:ceil((m2-m1)/m):m2,n1:ceil((n2-n1)/n):n2) + v(m1:ceil((m2-m1)/m):m2,n1:ceil((n2-n1)/n):n2) .* v(m1:ceil((m2-m1)/m):m2,n1:ceil((n2-n1)/n):n2)).^0.5);
        pcolor(X,Y,norm_u');
        shading interp;
        u_t = u(m1:ceil((m2-m1)/m):m2,n1:ceil((n2-n1)/n):n2) ./ norm_u;
        v_t = v(m1:ceil((m2-m1)/m):m2,n1:ceil((n2-n1)/n):n2) ./ norm_u;
        quiver(X,Y,u_t',v_t','black')
        
        str = sprintf('%i',K0);
        saveas(gcf,[fpat,filesep,str],'jpg');
        close all;
    end
    
    
    
    if(Daten)
        load('Matrizen.mat');        
    else
        disp('Init A');
        tic
        %% Divergenzmatrix aufstellen
        A = sparse((M-1)*(N-1),2*M*N);
        for i=1:M-1
            for j=1:N-1
                % Zelle i,j
                A((i-1)*(N-1)+j,(i-1)*(N)+j)              = -1;
                A((i-1)*(N-1)+j,(i-1)*(N)+j+1)            = -1;
                A((i-1)*(N-1)+j,(i  )*(N)+j)              =  1;
                A((i-1)*(N-1)+j,(i  )*(N)+j+1)            =  1;

                A((i-1)*(N-1)+j,(M)*(N)+(i-1)*(N)+j)      = -1;
                A((i-1)*(N-1)+j,(M)*(N)+(i  )*(N)+j)      = -1;
                A((i-1)*(N-1)+j,(M)*(N)+(i-1)*(N)+j+1)    =  1;
                A((i-1)*(N-1)+j,(M)*(N)+(i  )*(N)+j+1)    =  1;
            end
        end

        disp('Zeit zur Initialisierung von A')
        time = toc

        save MatrizenWeniger A L U P Q
    end
    

    %% Lösungsmethode vorbereiten
    A_NS = sparse(A(:,vNS));
    if(direkt)
        tic
        Z = (A_NS * A_NS' + alpha * speye((M-1)*(N-1),(M-1)*(N-1))); 
        disp('Zeit für Produkt');
        time = toc
        disp('Zeit für Sparsifizierung')
        tic
        Z = sparse(Z);
        time = toc
        disp('Zeit für Zerlegung')
        tic
        [L,U,P,Q] = lu(Z);
        %[R,P,p] = chol(Z);
        time = toc
    else
%        W = sparse( [eye(size(A_NS,1),size(A_NS,1)),A_NS;A_NS',zeros(size(A_NS,2),size(A_NS,2))] );
        P = sparse( A_NS*A_NS' + alpha*eye((M-1)*(N-1),(M-1)*(N-1)) );
    end    
    
    
    %% Hilfsvariablen
    lap_u       = zeros(size(u));
    lap_v       = zeros(size(v));
    x_u         = zeros(size(u));
    x_v         = zeros(size(v));
    x_phi_links = zeros(size(v));
    x_phi_rechts= zeros(size(v));
    y_u         = zeros(size(u));
    y_v         = zeros(size(v));
    y_phi_unten = zeros(size(v));
    y_phi_oben  = zeros(size(v));
    mit_u       = zeros(size(u));
    mit_v       = zeros(size(v));
    mit_phi     = zeros(size(v));
    c           = zeros(2*M*N,1);
    
    sum_time = 0;
    stream_on = true;
    stream_count = 0;
    %% Zeitschritte
    for k = 1:K
        stream_count = stream_count + 1;
        if(mod(stream_count,pt_load) == 0 && stream_on)
            stream_on = false;
            stream_counter = 0;
        elseif(mod(stream_count,pt_empt) == 0 && stream_on == false)
            stream_on = true;
            stream_counter = 0;
        end
        disp('-----------------------')
        tic
%         u_in    = u_in + 0.0001;
%         v_in    = v_in + 0.0001;
        
        str = sprintf('Schritt %i, dt = %d',k,dt);
        
        disp(str);
        
        
        %% Hilfsvariablen berechnen
        % Laplace-Operator
        lap_u(2:M-1,2:N-1)      = (0.25/(dx^2)) * ( u(1:M-2,2:N-1)+u(3:M,2:N-1)+u(2:M-1,1:N-2)+u(2:M-1,3:N)-4*u(2:M-1,2:N-1) );
        lap_v(2:M-1,2:N-1)      = (0.25/(dx^2)) * ( v(1:M-2,2:N-1)+v(3:M,2:N-1)+v(2:M-1,1:N-2)+v(2:M-1,3:N)-4*v(2:M-1,2:N-1) );
        
        % x-Ableitung
        x_u(2:M-1,2:N-1)         = (0.5/dx) * ( -u(1:M-2,2:N-1)+u(3:M,2:N-1) );
        x_v(2:M-1,2:N-1)         = (0.5/dx) * ( -v(1:M-2,2:N-1)+v(3:M,2:N-1) );
        
        % y-Ableitung
        y_u(2:M-1,2:N-1)         = (0.5/dx) * ( -u(2:M-1,1:N-2)+u(2:M-1,3:N) );
        y_v(2:M-1,2:N-1)         = (0.5/dx) * ( -v(2:M-1,1:N-2)+v(2:M-1,3:N) );
        
        % Mittel der vier Nachbarn
        mit_u(2:M-1,2:N-1)      = 0.25 * ( u(1:M-2,2:N-1)+u(3:M,2:N-1)+u(2:M-1,1:N-2)+u(2:M-1,3:N) );
        mit_v(2:M-1,2:N-1)      = 0.25 * ( v(1:M-2,2:N-1)+v(3:M,2:N-1)+v(2:M-1,1:N-2)+v(2:M-1,3:N) );
        mit_phi(2:M-1,2:N-1)      = 0.25 * ( phi(1:M-2,2:N-1)+phi(3:M,2:N-1)+phi(2:M-1,1:N-2)+phi(2:M-1,3:N) );
        
        %% Zeitschritt
        u = mit_u + dt * (- mit_u .* x_u - mit_v .* y_u + ny * lap_u);
        v = mit_v + dt * (- mit_u .* x_v - mit_v .* y_v + ny * lap_v);
        
        %% Randbedingungen
        % Einstrom-Flächen
        u(1,:)      = u_in * ones(1,N);
        %u(1,110:140)      = u_in * ones(1,31);
        u(M,:)      = u(M-1,:);
        v(1,:)      = 0 * ones(1,N);
        v(M,:)      = 0 * ones(1,N);
        % No-Slip-Ränder
        u(:,1)      = u_in * ones(M,1);
        u(:,N)      = u_in * ones(M,1);
        v(:,1)      = 0 * ones(M,1);
        v(:,N)      = 0 * ones(M,1);
        % Einbau
        u(Einbau)   = Null(Einbau);
        v(Einbau)   = Null(Einbau);
        
        %% Divergenzbereinigung
        %  Das Geschwindigkeitsfeld wird auf einen Schlag um seine
        %  Divergenz bereinigt.
        % Lösung von Matrix nach Vektor
        for i=1:M
            c((i-1)*N+1:i*N,1) = u(i,:)';
        end
        for i=1:M
            c(M*N+(i-1)*N+1:M*N+i*N,1) = v(i,:)';
        end
        
        % dt = 1e-3 * 1/max(c);
        
        % Residuum bestimmen
        r = A * c;
        % Divergenzkorrektur im Least-Squares-Sinne
        %dc = - A_NS' * ((A_NS * A_NS' + 0.1 * eye((M-1)*(N-1),(M-1)*(N-1)))\r);
        %dc = - A_NS' * (Z\r);
        if(direkt)
            dc = A_NS' * ( - Q * ( U\ ( L\ (P*r) ) ) );
        else
%             % Iterative Lösung über Symmetrisches Ersatzsystem
%             % A = sparse((M-1)*(N-1),2*M*N);
%             b = [-r;zeros(size(A_NS,2),1)];
%             h = minres(W,b,tol);
%             dc = A_NS' * h(size(r,1)+1:size(h,1),1);
            h = pcg(P,r,1e-2,100);
            dc = A_NS' * h;
        end
        % Update um LS-Minimalkorrektur
        c(vNS) = c(vNS) + dc;
        % Lösung von Vektor nach Matrix
        for i=1:M
            u(i,:) = c((i-1)*N+1:i*N,1)';
        end
        for i=1:M
            v(i,:) = c(M*N+(i-1)*N+1:M*N+i*N,1)';
        end
        
        %% Zeitschritt für Steam mittels Upwind
        phi(2:M-1,2:N-1) = mit_phi(2:M-1,2:N-1) - 0.5*dt * ( u(3:M,2:N-1) .* phi(3:M,2:N-1) - u(1:M-2,2:N-1) .* phi(1:M-2,2:N-1) + v(2:M-1,3:N) .* phi(2:M-1,3:N) - v(2:M-1,1:N-2) .* phi(2:M-1,1:N-2));
        
        
        %% Rauch
        % Rand
        phi(1,:)    = Null(1,:);
        phi(M,:)    = Null(M,:);
        phi(:,1)    = Null(:,1);
        phi(:,N)    = Null(:,N);
        if(stream_on)
            % Eingangsdüsen
            phi(Duese1)  = Eins(Duese1);
        end
        %Düse2
        phi(Duese2)  = Eins(Duese2);
        % Einbau
        phi(Einbau) = Null(Einbau);
        
        disp('Zeit für diesen Schritt');
        time = toc
        disp('geschätzte Restzeit')
        sum_time = sum_time + time;
        Restzeit = sum_time/k * (K-k)
        
        
        %% Zeitschrittfeld plotten
        if(mod(k,dK) == 0)
            str =  sprintf('Schritt %i, dt = %d',k+K0,dt);
%             figure;
%             hold on;
%             colormap((jet+white)/2);
%             grid off;
%             axis equal;
%             title(str);
%             xlabel('x');
%             ylabel('y');
%             norm_u = ((u .* u + v .* v).^0.5);
%             pcolor(norm_u');
%             shading interp;
%             u_t = u ./ norm_u;
%             v_t = v ./ norm_u;
%             quiver(u_t',v_t','black');
            
            %% Für Detail-Ansichten
            m1 = 1
            m2 = 400
            m  = 100
            n1 = 1
            n2 = 250
            n  = 62

            clf;
            colormap((jet+white)/2);
            hold on;
            grid off;
            axis equal;
            title(str);
            xlabel('x');
            ylabel('y');

            [X,Y] = meshgrid(m1:ceil((m2-m1)/m):m2, n1:ceil((n2-n1)/n):n2);
            
            norm_u = ((u(m1:ceil((m2-m1)/m):m2,n1:ceil((n2-n1)/n):n2) .* u(m1:ceil((m2-m1)/m):m2,n1:ceil((n2-n1)/n):n2) + v(m1:ceil((m2-m1)/m):m2,n1:ceil((n2-n1)/n):n2) .* v(m1:ceil((m2-m1)/m):m2,n1:ceil((n2-n1)/n):n2)).^0.5);
            GG = max(norm_u);
            GG = max(GG);
            dt = 0.9/GG;
%             pcolor(X,Y,norm_u');
            
            colormap bone;
            pcolor([1:M],[1:N],-phi');
            
            shading interp;
            u_t = u(m1:ceil((m2-m1)/m):m2,n1:ceil((n2-n1)/n):n2) ./ norm_u;
            v_t = v(m1:ceil((m2-m1)/m):m2,n1:ceil((n2-n1)/n):n2) ./ norm_u;
            quiver(X,Y,u_t',v_t','black')
            
            str = sprintf('%i',k+K0);
            saveas(gcf,[fpat,filesep,str],'jpg');
            
        end
    end
    
    if(false)
        %% Endzeitpunkt plotten
        figure;
        title('Endzeitpunkt');
        colormap((jet+white)/2);
        hold on;
        grid off;
        axis equal;
        title(str);
        xlabel('x');
        ylabel('y');
        pcolor(((u .* u + v .* v).^0.5)');
        shading interp;
        [X,Y] = meshgrid(1:ceil(M/30):M, 1:ceil(N/30):N);
        quiver(X,Y,u(1:ceil(M/30):M, 1:ceil(N/30):N)',v(1:ceil(M/30):M, 1:ceil(N/30):N)','black')
    end
    
    if(false)
        %% Für Detail-Ansichten
        m1 = 90
        m2 = 220
        m  = 100
        n1 = 120
        n2 = 190
        n  = 50

        figure;
        title('Endzeitpunkt');
        colormap((jet+white)/2);
        hold on;
        grid off;
        axis equal;
        title(str);
        xlabel('x');
        ylabel('y');

        [X,Y] = meshgrid(m1:ceil((m2-m1)/m):m2, n1:ceil((n2-n1)/n):n2);
        norm_u = ((u(m1:ceil((m2-m1)/m):m2,n1:ceil((n2-n1)/n):n2) .* u(m1:ceil((m2-m1)/m):m2,n1:ceil((n2-n1)/n):n2) + v(m1:ceil((m2-m1)/m):m2,n1:ceil((n2-n1)/n):n2) .* v(m1:ceil((m2-m1)/m):m2,n1:ceil((n2-n1)/n):n2)).^0.5);
        pcolor(X,Y,norm_u');
        shading interp;
        u_t = u(m1:ceil((m2-m1)/m):m2,n1:ceil((n2-n1)/n):n2) ./ norm_u;
        v_t = v(m1:ceil((m2-m1)/m):m2,n1:ceil((n2-n1)/n):n2) ./ norm_u;
        quiver(X,Y,u_t',v_t','black')
    end
    
end