function r = lindberg_resid ( t, y, dydt )

%*****************************************************************************80
%
%% lindberg_resid() returns the residual of the Lindberg ODE.
%
%  Discussion:
%
%    Note that components y1(t) and y2(t) first sink to extraordinarily
%    small values, and then undergo explosive growth, sometime before t=1.
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    11 May 2024
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Bengt Lindberg,
%    On a dangerous property of methods for stiff differential equations,
%    BIT Numerical Mathematics,
%    Volume 14, 1974, pages 430-436.
%
%    Daniel Watanabe, Qasim Sheikh,
%    One-leg formulas for stiff ordinary differential equations,
%    SIAM Journal on Scientific and Statistical Computing,
%    Volume 5, Number 2, June 1984, pages 489-496.
%
%  Input:
%
%    real T(N): selected time values.
%
%    real Y(N,4): corresponding solution values.
%
%    real DYDT(N,4): corresponding derivate values.
%
%  Output:
%
%    real R(N,4): the residuals.
%
  n = length ( t );
  r = zeros ( n, 4 );

  r(:,1) = dydt(:,1) - (   10.0^4 * y(:,1) .* y(:,3) + 10.0^4 * y(:,2) .* y(:,4) );  
  r(:,2) = dydt(:,2) - ( - 10.0^4 * y(:,1) .* y(:,4) + 10.0^4 * y(:,2) .* y(:,3) ); 
  r(:,3) = dydt(:,3) - (   1.0 - y(:,3) );
  r(:,4) = dydt(:,4) - ( - 0.5 * y(:,3) - y(:,4) + 0.5 );
 
  return
end

