function [x,y]=back_euler(f,xRange,yInitial,numSteps)
% [x,y]=back_euler(f,xRange,yInitial,numSteps) computes
% the solution to an ODE by the backward Euler method
%
% f is a function handle for a function with signature
%   [fValue fPartial]=f(x,y)
%   where fValue is a column vector and fPartial a matrix
% xRange is a two dimensional vector of beginning and
%   final values for x
% yInitial is a column vector for the initial value of y
% numSteps is the number of evenly-spaced steps to divide
%   up the interval xRange
% x is a row vector of selected values for the 
%   independent variable
% y is a matrix. The k-th column of y is 
%   the approximate solution at x(k)

% M. M. Sussman
% $Id: back_euler.m,v 1.4 2017/01/16 01:14:20 mike Exp $

% force x to be a row vector
x(1,1) = xRange(1);
h = ( xRange(2) - xRange(1) ) / numSteps;
y(:,1) = yInitial;
for k = 1 : numSteps
  x(1,k+1) = x(1,k) + h;

  Y = (y(:,k)) + h * feval( f, x(1,k), y(:,k));
  [Y,isConverged]= newton4euler(f,x(k+1),y(:,k),Y,h);

  if ~ isConverged
    error(['back_euler failed to converge at step ', ...
          num2str(k)])
  end

  y(:,k+1) = Y;
end
