function [J] = jacobian(F,t,Y,h)
  
% numerical jacobian of g(y(i+1)) = y(i+1) - y(i) - (hF(t(i+1),y(i+1)))
%input: 
%       F = @(x) function handle 
%       t = a single value
%       Y = column vector 
%       h = stepsize  

  Y = Y';
  eps = 1e-5;
  X = [t Y];
  X_d = X;
  fx = feval(F,t,Y);
  n = length(X);
  for j=1:n
    X_d(j) = X_d(j) + eps;
    J(:,j) = (feval(F,X_d(1),X_d(2:n))-fx)/eps;
    X_d(j) = X(j);
    end
    J = J(:,2:n);
    n = size(J,1);
    I = eye(n);
    J = I - h.*J;
    
    
  end