function [x, y] = ImpEuler( x0, x_end, y0, h, rhs_fctn )%%  Purpose:%  Solve the initial value problem%               %                      y'(x) = f(x,y(x))%                      y(x0) = y0%%  on the interval [ x0, x_end ] using the implicit Euler method%  The implicit Euler method is given by%%    y(:,i+1) - h * f(x(i+1),y(:,i+1)) - y(:,i) = 0%%%%  Usage:%%       [x, y] = ImpEuler( x0, x_end, y0, h, rhs_fct, rhs_fct_y )%%  Parameters:%%  Input:%%  x0          initial x%%  x_end       final x%%  y0          initial values%%  h           step size%%  rhs_fctn    function evaluating f(x,y) for given x and y%              rhs_fctn(x, y, '') returns  f(x,y)%              rhs_fctn(x, y, 'jacobian') returns Jacobian df(x,y)/dy%%%  Output:%  %  x        vector of x values x(i) = x0 + (i-1)*h%%  y        approximation of solution at x%% find number of time steps needed to reach final timemx = ceil((x_end-x0)/ h);% allocate workspacen = size(y0(:),1);x = (x0:h:x0+mx*h);y = zeros(n,mx+1);% tolerance for Newton's methodNewt_tol = 0.01*h^2;% Set initial valuesy(:,1) = y0(:);for i = 1:mx    % solve the nonlinear system for y(:,i+1) using Newton's     % method with starting value  y(:,i)    y(:,i+1) = y(:,i);  % starting value for Newton's method    Newt_it  = 0;       % counter for Newton iterations    Newt_fct = y(:,i+1) - h * feval(rhs_fctn,x(i+1),y(:,i+1),'') - y(:,i);    while( norm(Newt_fct) > Newt_tol*norm(y(:,i+1)) )        % Compute Jacobian         Newt_Jac  = eye(n) - h * feval(rhs_fctn,x(i+1),y(:,i+1),'jacobian');        Newt_step = - (Newt_Jac \ Newt_fct);   % Newton step        y(:,i+1)  = y(:,i+1) + Newt_step;      % new Newton iterate        Newt_it   = Newt_it + 1;              Newt_fct  = y(:,i+1) - h * feval(rhs_fctn,x(i+1),y(:,i+1),'') - y(:,i);    end    %disp([ x(i+1), Newt_it ])end