function [u,xc] = entrausch_nl(n,pn,mu)

% das Gitter
    h  = 2/n;                                      % Zellenbreite
    xc = -1.0 + 0.5*h + h*(0:(n-1))';              % Zellzentren
    xg = -1.0         + h*(0:n)';                  % Zellschnittstellen

% verrauschte Funktion
    k  = fix(n/2);
    v0 = [zeros(k,1);ones(n-k,1)];                 % Treppenfunktion

%   v0 = sin(10*xc);                               % mehrere Treppen
%   v0 = double(v0 > 0) - double(v0 <= 0);

    randn('seed',0);                               % init randn
    v  = v0 + pn*randn(n,1);                       % verrauschte Funktion

% Begin mit linearer Abschaetzung
    D  = diag([ones(n-1,1);0])...                  % Diskretisierung
       + diag([0;ones(n-1,1)])...                  %   der zweiten
       - diag(ones(n-1,1),1)...                    %   Ableitung mit
       - diag(ones(n-1,1),-1);                     %   Randbedingungen
    A  = eye(n) + mu*D/h^2;                        % entrauschender Operator
    u  = A\v;                                      % entrauschte Funktion
    plot(xc,v0,'b',xc,v,'r',xc,u,'g');             % grafisch darstellen

% Picard Iteration: setzte fort bis Unterschiede klein genug
    itmax = 100;                                   % max Iterationen
    tol   = 1.0e-4;                                % Abbruchkriterium
    omega = 1.0;                                   % Daempfung
    it    = 0;
    us    = v;
    while ((max(abs(u-us)) > tol*max(abs(u))) & (it < itmax))
        us = u;                                    % speichere altes u
        ux = (u(2:n)-u(1:(n-1)))/h;                % ux in Zellschnittstellen
        ps = phistrich(ux.^2);                     % Leitwerte
        D  = diag([ps;0])...                       % Diskretisierung des
           + diag([0;ps])...                       %   nicht linearen
           - diag(ps,+1)...                        %   entrauschenden
           - diag(ps,-1);                          %   Operators mit
        A  = eye(n) + mu*D/h^2;                    %   Randbedingungen
        u  = A\v;                                  % aktuelle Abschaetzung
        u  = us + omega*(u - us);                  % Daempfung
        it = it+1;
        plot(xc,v0,'b',xc,v,'r',xc,u,'g');         % grafisch darstellen
        drawnow;
    end

% zeige exakt, verrauscht und entrauscht zusammen
    plot(xc,v0,'b',xc,v,'r',xc,u,'g');             % grafisch darstellen

% Leitwertsfunktion:  Energieterm ist phi(ux^2).
function ps = phistrich(ux2)
    q  = 1;                                        % phi(s)=[s+ep]^(q/2)
    ep = 1.0e-2;                                   % phi'(s)=
    ps = (q/2)./(ux2 + ep).^(1-q/2);               %   (q/2)[s+ep]^(q/2-1)
