function J=dG_BD2(dt,U,Yn,Yp,Iinj)

lY=length(Yn);

deltaU=.1e-3;

for j=1:lY
    Up=U;
    Up(j)=Up(j)+deltaU;
    Fp=G_BD2(dt,Up,Yn,Yp,Iinj);
    J(:,j)=(Fp-G_BD2(dt,U,Yn,Yp,Iinj))/deltaU;
end


     
