procedure findu(var v: vector; var u: ltuple);
var
  i,j,k: integer;
  s1,s2,sk: scalar;
  fu,u0,du,lt,w: ltuple;
  dfu,m,mt: matrix;
begin { findu }
  write('(findu'); flush(output);
  for i:=0 to 3 do shrink(u[i]);

  for k:=1 to 10 do begin    { Newton for solving fu=0 approx }
    u0:=u;
    df(v,u,fu,dfu);
    minverse(dfu,m);
    for i:=0 to 3 do for j:=0 to 3 do shrink(m[i,j]);
    lmprod(m,fu,du);
    ldiff(u,du,lt);
    u:=lt;
    for i:=0 to 3 do shrink(u[i]);
  end;

  if normalize then begin
    for i:=0 to lmax do begin              { weights for norm }
      sabs(du[i],w[i]);
      w[i].l:=w[i].u;
    end;

    lball(stwo,w,lt);
    lsum(u0,lt,u);               { ball of radius 2 around u0 }
    df(v,u,fu,dfu);
    mprod(m,dfu,mt);
    for i:=0 to 3 do begin
      sdiff(mt[i,i],sone,s1);
      mt[i,i]:=s1;
    end;
    mnorm(w,mt,sk);              { norm of derivative on ball }
    writeln;

    if not sle(sk,shalf) then writeln('findu: error');
    sdiff(sone,sk,s1);
    squot(sone,s1,s2);
    lball(s2,w,lt);
    lsum(u0,lt,u);
  end;
  write(')'); flush(output);
end { findu };

