function J=jacob(FH,dt,U,Yn,Iinj)

lY=length(Y);

deltaY=.1e-5;

for j=1:lY
    Yp=Y;
    Yp(j)=Yp(j)+deltaY;
    Fp=FH(Yp,Iinj);
    J(:,j)=(Fp-FH(Y,Iinj))/deltaY;
end


     
