function [ t, u, v ] = pendulum_euler_backward ( n )

%*****************************************************************************80
%
%% pendulum_euler_backward uses the backward Euler method on the pendulum ODE.
%
%  Discussion:
%
%    We write down the backward Euler method for the pendulum:
%
%      u(i+1) = u(i) + dt *       v(i+1);
%      v(i+1) = v(i) - dt * g/l * u(i+1);
%
%    Bring the future terms to the left hand side.
%
%      u(i+1) - dt       * v(i+1) = u(i)
%      v(i+1) + dt * g/l * u(i+1) = v(i)
%
%    Rewrite as a matrix equation:
%
%      ( 1,      -dt ) * ( u(i+1) ) = u(i)
%      ( dt*g/l,  1  )   ( v(i+1) ) = v(i)
%
%    The matrix A is a kind of unnormalized rotation.
%
%      A = (  1,     -dt )
%          ( dt*g/l,   1 )
%
%    The inverse is:
%
%      inv(A) = 1/(1+g/l*dt^2) * ( 1,     dt )
%                                (-g/l*dt, 1 )
%
%    Multiply the matrix equation A*(u(i+1),v(i+1)) = (u(i),v(i))
%    by the inverse:
%
%      u(i+1) = 1.0 / ( 1.0 + g/l*dt^2 ) * (            u(i) + dt * v(i) );
%      v(i+1) = 1.0 / ( 1.0 + g/l*dt^2 ) * ( - g/l*dt * u(i) +      v(i) );
%
%    This is an explicit equation for the new values.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    27 February 2020
%
%  Author:
%
%    John Burkardt
%
%  Input:
%
%    integer N: the number of steps to take.
%
%  Output:
%
%    real t(n+1), u(n+1), v(n+1): the times and solution values.
%
  g = 9.8;
  l = 1.0;

  t = zeros ( n + 1, 1 );
  u = zeros ( n + 1, 1 );
  v = zeros ( n + 1, 1 );

  a = 0.0;
  b = 10.0;
  dt = ( b - a ) / n;

  t(1) = a;
  u(1) = 1.0;
  v(1) = 0.0;

  for i = 1 : n
    t(i+1) = t(i) + dt;
    u(i+1) = 1.0 / ( 1.0 + g/l * dt^2 ) * (            u(i) + dt * v(i) );
    v(i+1) = 1.0 / ( 1.0 + g/l * dt^2 ) * ( - g/l*dt * u(i) +      v(i) );
  end

  return
end

