procedure df(var v: vector; var u,fu: ltuple; var dfu: matrix);
var
  i: integer;
  v0,vt,vt1,vt2: vector;
begin { df }
  pexp(3,u[1],v0);                         { v0 ~= exp(u[1].) }
  vtrans(4,rho,u[3],v,vt);                  { vt ~= v(u[3]+.) }

  pscale(3,u[2],vt,vt1);
  pprod(3,v0,vt1,vt2);
  for i:=0 to 3 do begin
    sprod(u[0],vt2[i],vt1[i]);
    sdiff(vt1[i],lnormal[i],fu[i]);    
    dfu[i,0]:=vt2[i];
    if i=0 then szero(dfu[i,1]) else dfu[i,1]:=vt1[i-1];
  end;

  pder(3,vt,vt2);
  pscale(3,u[2],vt2,vt1);
  pprod(3,v0,vt1,vt2);
  for i:=0 to 3 do begin
    sprod(u[0],vt2[i],vt1[i]);
    if i=0 then szero(dfu[i,2]) else dfu[i,2]:=vt1[i-1];
    dfu[i,3]:=vt1[i];
  end;
end { df };

