function [yy, tt] = ode_bdf4_var (f, int, y0, tol, Nmax)

  t0   = int(1);
  t1   = int(2);
  h    = (t1-t0)/Nmax;      % initial step size
  hmax = (t1-t0)/10;
  hmin = (t1-t0)/1e10;
  y    = y0(:);             % Make sure we initialze with a column vector
  yy   = y;
  t    = t0;
  tt   = t0;
  step = 0;                 % Step counter
  m = 4;                    % Order of BDF method
  n = length(y0);           % Dimension of system
  
  % Initialize the Nordsieck representation
  %
  % We are performing a simple initialization of the code, where all
  % derivatives are initially to be taken as zero.  This corresponds to
  % using a first order method for the first step, which is OK since
  % step size control is going to keep the error small.

  % Note that each vector of Taylor coefficients is a ROW vector
  z = zeros (n,m+1);
  z(:,1) = y;

  % Nordsieck vector for order 4 (Coefficients computed on Mathematica) 
  v = [12/25 1 7/10 1/5 1/50];

  % Shift matrix S
  % This is the transpose of the usual definition, because we are
  % keeping Taylor coefficients as row vectors.
  S = zeros(m+1);
  for j = 0:m
    for k = 0:j
      S(1+j, 1+k) = bincoeff(j,k);
    end
  end

  % Unit vector e1
  e1 = zeros (m+1,1);
  e1(2) = 1;
  
  B = eye(n);          % approximate Jacobian
    
  while t(length(t)) < t1

    step = step+1;
    if step >= Nmax
       disp ('Maximum Number of steps reached before completion.');
       return;
    end

    if h < hmin
       disp ('The step size has become too small.  Singularity?');
       return;
    end
    
    yi = z(:,1);       % Old value of y 
    b  = z*S;          % Push Taylor coefficients forward
    ye = b(:,1);       % New value of y via extrapolation

    % Solve implicit step via Broyden iteration

    b  = b - (b * e1) * v;
    fi = b(:,1) + 12/25 * h * feval (f, t+h, yi) - yi;

    for k = 1:20
      dy = -B\fi;
      yi = yi + dy;
      if norm(dy) < 10e-11
         break;
      end
      fj = b(:,1) + 12/25 * h * feval (f, t+h, yi) - yi;
      df = fj - fi;
      fi = fj;
      B = B + ((df - B*dy) * dy')/(dy' * dy);     
    end

    % Check error and compute new stepsize

    error = norm (ye-yi);

    if error < tol
       % Accept Step
       z = b + h * feval (f, t+h, yi) * v;
       yy = [yy z(:,1)];
       t  = t+h;
       tt = [tt t];
       hn = min ([hmax, 2*h, 0.85*h*(tol/error)^(1/6)]);   
       if t(length(t)) + hn > t1
          hn = t1 - t(length(t));
       end
    else % Reject Step
       hn = h/2; 
    end

    % Adapt vector of Taylor coefficients to new step size
    z = z * diag((hn/h).^[0:m]);
    h = hn;
    
  end % Main loop
    
end