%
% 2 x 2 ODE system by variable order BDF
%
% Global variables
  global a b ncall;
%
% Model parameters
  a=500000.5;
  b=499999.5; 
%
% Select method
  for mf=1:2
%
% Error tolerances
  reltol=1.0e-02;
  abstol=1.0e-02;
%
% Cases for changes in error tolerances
  for ncase=1:2
    reltol=1.0e-02*reltol;
    abstol=1.0e-02*abstol;
%
%   Variables for ODE integration
    t0=0.0;
    tf=10.0;
    tout=[t0:1.0:tf]'; 
    nout=11;
%
%   Initialize number of calls to derivative
%   subroutine
    ncall=0;
%
%   Initial condition
    y10=0.0;
    y20=2.0;
    y0=[y10 y20]';
%
%   Call ODE integrator (caution: calling ode23, ode45 in place
%   of ode23s, ode15s in the following statements will result
%   in long computer runs for a=500000.5, b=499999.5; see the 
%   discussion in Appendix C for details)
    options=odeset('RelTol',reltol,'AbsTol',abstol);
    if(mf==1)[t,y]=ode23s('ode2x2',tout,y0,options); end    
    if(mf==2)[t,y]=ode15s('ode2x2',tout,y0,options); end 
%       
%   Display solution and error
    fprintf('\n\n  mf = %1d\n  case = %1d\n  reltol = %6.2e\n  abstol = %6.2e\n\n',...
            mf,ncase,reltol,abstol); 
    fprintf('    t      y1e       y1          erry1\n           y2e       y2          erry2\n');
    for i=1:nout
      lambda1=-(a-b);
      lambda2=-(a+b);
      exp1=exp(lambda1*t(i));
      exp2=exp(lambda2*t(i));
      y1e=(y10+y20)/2.0*exp1-(y20-y10)/2.0*exp2;
      y2e=(y10+y20)/2.0*exp1+(y20-y10)/2.0*exp2;
      erry1=y1e-y(i,1);
      erry2=y2e-y(i,2);
      fprintf('%5.1f%9.4f%9.4f%15.10f\n     %9.4f%9.4f%15.10f\n\n',...
              t(i),y1e,y(i,1),erry1,y2e,y(i,2),erry2);
    end
%
% Next case
  fprintf('  ncall = %5d\n',ncall);
  end
%
% Next method
  end
%
% Plot last solution
  plot(t,y);
  xlabel('t')
  ylabel('y1(t),y2(t)')
  title(' Appendix C, 2 x 2 Linear System')
  gtext('y1(t)');
  gtext('y2(t)');
  print appc.ps


