function [u,x] = regis(n,pn,mu)

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

% Funktionen zu registrieren
    k1 = fix(n/8);
    k2 = fix(3*n/8);
    I0 = sin((1:(n-k1-k2))/10)';
    I0 = [zeros(k1,1);I0;zeros(k2,1)];             % 1. Funktion
    k1 = fix(3*n/8);
    k2 = fix(n/8);
    I1 = sin((1:(n-k1-k2))/10)';
    I1 = [zeros(k1,1);I1;zeros(k2,1)];             % 2. Funktion

% verrauschte Funktionen
    randn('seed',0);                               % init randn
    I0 = I0 + pn*randn(n,1);                       % 1. verrauschte Funktion
    I1 = I1 + pn*randn(n,1);                       % 2. verrauschte Funktion

% Begin mit Identitaet
    u  = x;

% Diskretisierung der vierten Ableitung mit Randbedingungen
    D = 6*diag(ones(n  ,1)   ) ...
      - 4*diag(ones(n-1,1),-1) ...
      - 4*diag(ones(n-1,1),+1) ...
      +   diag(ones(n-2,1),-2) ...
      +   diag(ones(n-2,1),+2);

    D(  1,  1) =  1;
    D(  n,  n) =  1;

    D(  1,  2) = -2;
    D(  2,  1) = -2;
    D(  n,n-1) = -2;
    D(n-1,  n) = -2;

    D(  2,  2) =  5;
    D(n-1,n-1) =  5;

% Ableitung von I1 in x
    Ix = [2*I1(1)-I1(2);I1;2*I1(n)-I1(n-1)];
    Ix = (Ix(3:(n+2))-Ix(1:n))/h;

% Quasi-Newton Iteration:  setze fort bis Unterschiede klein genug
    itmax = 1000;                                  % max Iterationen
    tol   = 1.0e-4;                                % Abbruchkriterium
    omega = 1.0;                                   % Daempfung
    it    = 0;
    us    = u+1;
    while ((max(abs(u-us)) > tol*max(abs(u))) & (it < itmax))
        us = u;                                    % speichere altes u
        Is = interp1(x,Ix,u,'linear','extrap');    % Ableitung von I1 in u
        Iu = interp1(x,I1,u,'linear','extrap');    % I1 in u
        dI = (I0 - Iu).*Is;                        % treibende Kraft
        v  = -(mu/h^4)*D*u + dI;                   % Residuum
        A  = (mu/h^4)*D + diag(Is.^2);             % Residuum-Ableitung
        du = A\v;                                  % Newton Schritt
        u  = u + omega*du;                         % gedaempfte Abschaetzung
        it = it+1;
        plot(x,I0,'b',x,I1,'g',x,Iu,'r');          % grafisch darstellen
        drawnow;
        pause(0.01);
    end

% interpoliertes Ergebnis grafisch darstellen
    nb = n;
    for i=1:nb
        z  = (i-1)/(nb-1);                         % interpolierte Funktion
        Iz = interp1(x,I1,x+z*(u-x),'linear','extrap')';
        plot(x,I0,'b',x,I1,'g',x,Iz,'r');          % grafisch darstellen
        pause(0.01);
    end
