%
% 2 x 2 ODE system by first order BDF
%
% Open a file for output
  fid=fopen('appc1.out','w');
%
% Number of ODEs
  neqn=2;
%
% Convergence tolerance, maximum number of
% iterations to compute Newton corrections
  eps=0.00001;
  maxiter=20; 
%
% Initial condition
  y10=0.0;
  y20=2.0;
  y(1)=y10;
  y(2)=y20;
  t=0.0;
%
% Values just for initial output
  dy(1)=0.0;
  dy(2)=0.0;
  niter=0;
% 
% Integration step, steps/output, number of outputs
  h=0.01;
  nsteps=100;
  nout=11;
%
% Problem parameters
  a=500000.5;
  b=499999.5;
%
% Jacobian matrix
  J(1,1)=-a*h-1.0;
  J(1,2)=b*h;
  J(2,1)=b*h;
  J(2,2)=-a*h-1.0;
%
% Print a heading
  fprintf(fid,'       t       dy(1)     dy(2)      y(1)      y(2)   erry(1)   erry(2)  iter\n');
  fprintf(    '       t       dy(1)     dy(2)      y(1)      y(2)   erry(1)   erry(2)  iter\n');
%
% nout outputs
  for i=1:nout
%
%   Initial output or after nstep integration steps completed;
%   display Newton corrections, numerical and exact solutions for
%   output
    lambda1=-(a-b);
    lambda2=-(a+b);
    exp1=exp(lambda1*t);
    exp2=exp(lambda2*t);
    y1e=exp1-exp2;
    y2e=exp1+exp2;
    erry1=y1e-y(1);
    erry2=y2e-y(2);
    fprintf(fid,'%10.5f%10.5f%10.5f%10.5f%10.5f%10.5f%10.5f%5d\n',...
                t,dy(1),dy(2),y(1),y(2),erry1,erry2,niter);
    fprintf(    '%10.5f%10.5f%10.5f%10.5f%10.5f%10.5f%10.5f%5d\n',...
                t,dy(1),dy(2),y(1),y(2),erry1,erry2,niter);
%
%   nsteps steps/output
    for is=1:nsteps
%
%     Initialize iteration counter
      niter=1;
      stop=0;
%
%     Test for the end of the current step
      while stop==0
%   
%       Functions g1, g2
        g=[h*(-a*y(1)+b*y(2))-y(1)+y10
           h*( b*y(1)-a*y(2))-y(2)+y20];
%
%       Solve for Newton corrections
%
%       Gaussian elimination
        dy=-J\g;
%
%       Jacobian inverse
%       dy=-J^(-1)*g;
%
%       Update solution
        for ic=1:neqn
          y(ic)=y(ic)+dy(ic);
        end   
        stop=1;     
%
%       Check if the corrections are within the tolerance eps
        for ic=1:neqn
          if abs(dy(ic))>eps
%
%           Convergence not achieved; continue calculation
            niter=niter+1;
            stop=0;
            break;
          end 
        end 
%     
%       If maximum iterations reached, accept current step
        if (niter==maxiter) stop=1; end
%
%     Continue integration step
      end
%
%     Integration step completed
      y10=y(1);
      y20=y(2);
      t=t+h;
%
%     Continue nstep integration steps
      end
%
% Continue integration for next output interval
  end
