function df=pdehypdf(t,u)
%PDEHYPDF Jacobian for HYPERBOLIC.

%       Copyright 1994-2014 The MathWorks, Inc.

  global pdehyp
  
  if ~(pdehyp.vq || pdehyp.vg || pdehyp.vh || pdehyp.vr || pdehyp.vc || pdehyp.va || pdehyp.vf)
    df=-pdehyp.K;
    return
  end
  
  nu=pdehyp.numConstrainedEqns;
  if(nu == pdehyp.totalNumEqns)
    uFull = u(1:nu);
  else
    uFull = pdehyp.B*u(1:nu) + pdehyp.ud;
  end
  
  if(pdehyp.vq || pdehyp.vh)
    [Q,G,H,R]=assemb(pdehyp.b,pdehyp.p,pdehyp.e,uFull,t);
  end
  
  if(pdehyp.vc || pdehyp.va)
    [K,M]=pdehyp.thePde.assema(pdehyp.p,pdehyp.t,pdehyp.c,pdehyp.a,zeros(pdehyp.N,1),uFull,t);
  end
  
  if pdehyp.vq
    pdehyp.Q=Q;
  end
  
  if pdehyp.vh
    pdehyp.B=pdenullorth(H);
  end
  
  if pdehyp.vc
    pdehyp.K=K;
  end
  
  if pdehyp.va
    pdehyp.M=M;
  end
  
  nu=size(pdehyp.B,2);
  if ~pdehyp.vh
    df=-[sparse(nu,nu) -speye(nu,nu); ...
      pdehyp.B'*(pdehyp.K+pdehyp.M+pdehyp.Q)*pdehyp.B sparse(nu,nu)];
  else
    dt=pdehyp.ts*eps^(1/3);
    t1=t+dt;
    dt=t1-t;
    tm1=t-dt;
    if pdehyp.vd
      [~,pdehyp.MM]=pdehyp.thePde.assema(pdehyp.p,pdehyp.t,0,pdehyp.d,zeros(pdehyp.N,1),uFull,t);
    end
    [Q,G,H,R]=assemb(pdehyp.b,pdehyp.p,pdehyp.e,uFull,t1);
    N1=pdenullorth(H);
    [Q,G,H,R]=assemb(pdehyp.b,pdehyp.p,pdehyp.e,uFull,tm1);
    Nm1=pdenullorth(H);
    BDot = (N1-Nm1)/(2*dt);
    BDotDot = (N1-2*pdehyp.B+Nm1)/(dt^2);
    if(isempty(pdehyp.C))
      nrFull = size(uFull,1);
      dampMat = sparse(nrFull, nrFull);
    else
      dampMat = pdehyp.C;
    end
    df=-[sparse(nu,nu) -speye(nu,nu); ...
      pdehyp.B'*((pdehyp.K+pdehyp.M+pdehyp.Q)*pdehyp.B+ ...
      pdehyp.MM*BDotDot+dampMat*BDot) ...
      2*pdehyp.B'*pdehyp.MM*BDot];
  end
  
end

