function df=pdeprbdf(t,u)
%PDEPRBDF Jacobian for PARABOLIC.

%       Copyright 1994-2014 The MathWorks, Inc.

global pdeprb

if ~(pdeprb.vq || pdeprb.vg || pdeprb.vh || pdeprb.vr || pdeprb.vc || pdeprb.va || pdeprb.vf)
  df=-pdeprb.K;
  return
end

if(pdeprb.numConstrainedEqns == pdeprb.totalNumEqns)
  uFull = u;
else
  uFull = pdeprb.B*u + pdeprb.ud;
end

if(pdeprb.vq || pdeprb.vh)
  [Q,G,H,R]=assemb(pdeprb.b,pdeprb.p,pdeprb.e,uFull,t);
end

if(pdeprb.vc || pdeprb.va)
  [K,M]=pdeprb.thePde.assema(pdeprb.p,pdeprb.t,pdeprb.c,pdeprb.a,zeros(pdeprb.N,1),uFull,t);
end

if pdeprb.vq
  pdeprb.Q=Q;
end

if pdeprb.vh
  pdeprb.B=pdenullorth(H);
end

if pdeprb.vc
  pdeprb.K=K;
end

if pdeprb.va
  pdeprb.M=M;
end

if ~pdeprb.vh
  df=-pdeprb.B'*(pdeprb.K+pdeprb.M+pdeprb.Q)*pdeprb.B;
else
  dt=pdeprb.ts*sqrt(eps);
  t1=t+dt;
  dt=t1-t;
  if pdeprb.vd
    [~,pdeprb.MM]=pdeprb.thePde.assema(pdeprb.p,pdeprb.t,0,pdeprb.d,zeros(pdeprb.N,1),uFull,t);
  end
  [Q,G,H,R]=assemb(pdeprb.b,pdeprb.p,pdeprb.e,uFull,t1);
  N1=pdenullorth(H);
  df=-pdeprb.B'*((pdeprb.K+pdeprb.M+pdeprb.Q)*pdeprb.B + ...
    pdeprb.MM*(N1-pdeprb.B)/dt);
end

end

