function [x,y] = BkEulerAuto(dfun,xspan,y0,AbsTol,RelTol)
% [x,y] = BkEulerFixH('dfun',tspan,y0,AbsTol,RelTol)
%
% Solves stiff systems of differential equations y' = f(x,y)
% by the backward Euler method with automatic step selection.
% 
% Arguments:
%   dfun    String, name of a function M-file
%             feval(dfun,x,y,'')         returns derivative
%             feval(dfun,x,y,'jacobian') returns jacobian
%   xspan   Integration interval endpoints [x0 xf] 
%   y0      Initial conditions at x0
%   AbsTol,RelTol [Optional] Error tolerances

% Set default tolerances if not specified by user.
% 
if nargin < 5
  RelTol = 1e-3;
  if nargin < 4
    AbsTol = 1e-6;
  end
end
x0 = xspan(1);                 % Initial point
xf = xspan(2);                 % Final point
dr = sign(xf-x0);              % Direction of integration
y0 = y0(:);                    % Initial values -> Column vector
neq = length(y0);              % # of dep. vars.
x = [x0; zeros(32,1)];         % Allocate space for steps
y = [y0'; zeros(32,neq)];      % Allocate space for solution
n  = 1;                        % Points to last step completed
init = 1;                      % boolean: starting integration
cvfail = 0;                    % boolean: Newton convergence failed
rtol = 1e-6;                   % Tolerance for Newton iteration
f0  = feval(dfun,x0,y0,'');    % derivative at initial point

% MAIN LOOP: Take steps until end of interval reached.

while dr*(xf-x) > 0
  if init

%   First Step: Estimate 2nd deriv bound with
%    steps 1/8, 1/64, 1/512 interval

    halfM2 = sqrt(eps);
    for p = 1:3
      h = (xf-x0)/8^p;
      y1 = y0 + h*f0;
      f1 = feval(dfun,x0+h,y1);
      E = norm( (h/2)*(f1-f0), inf);
      halfM2 = max([halfM2 E/h^2]);
    end

%   Select initial stepsize -- choose a step that
%   is expected to pass the error test, but is at
%   most 1/8 of the integration interval.

    h = sqrt(0.5*(AbsTol+RelTol*norm(y1,inf))/halfM2);
    h = min([h abs(xf-x0)/8]);
    h = dr*h;
    hfn = h*f0;
    init = 0;
  else

%   Not first step.
%   Select step based on error estimate in step just completed.

    if cvfail
      up = 0.25;
    else
      up = sqrt(0.5*(AbsTol+RelTol*norm(y(n,:),inf))/Est);
    end
    h = up*h;       % Scale the step
    hfn = up*hfn;   % ... and the saved derivative
  end

% Predict solution by forward Euler 
% [except on first step!] evaluate residual

  xnp = x(n) + h;
  if n>1
    z = hfn;
  else
    z = zeros(neq,1);
  end

% Evaluate and factor the iteration matrix for this step

  J = feval(dfun,x(n),y(n,:)','jacobian');
  [L,U] = lu(eye(neq)-h*J);

% Simplified Newton iteration.

  r = 1+rtol; m = 0; mxit = 10; cvfail = 0;
  while (norm(r) > rtol & m < mxit)
    r  = z - h*feval(dfun,xnp,y(n,:)'+z,'');
    dz = U \ ( L \ r );
    z  = z - dz;
    m = m+1;
  end
  if norm(r) <= rtol

%   Converged. Estimate error.

    Est = 0.5*norm(z-hfn,inf);
    if Est < AbsTol+RelTol*norm(y(n,:),inf)

%     Passed error test -- accept step, advance solution

      n = n+1;
      x(n)   = xnp;
      y(n,:) = y(n-1)+z';
      hfn = z;
      if mod(n,32)==0     % Allocate more space for solution
        x = [x; zeros(32,1)];
        y = [y; zeros(32,neq)];
      end
    end
  else
    
%   Newton convergence failed.
    cvfail = 1;
  end
end
x = x(1:n);      % Chop off unused array space
y = y(1:n,:);