% Script to compare erfc(x) with its asymptotic series for large x
%  erfc(x) = 2/sqrt(pi) int_x^infty exp(-y^2) dy
%          = exp(-x^2)/sqrt(pi) *(1/(2*x) - 1/(2^2*x^3) + 1*3/(2^3* x^5)...
% This can be written as the large-P limit of partial sums
%  s_P(x)  = sum_n=0^P t_n(x)
%  t_n(x)  = 2*exp(-x^2)/sqrt(pi) * a_n x^-(2*n+1)
% where
%  a_n = -a_{n-1} * (n-1/2), a_0 = 1/2.

  % Input parameters
  
  x = 0.5:0.1:4;
  x1 = 2;
  x2 = 3;
  P=10;

  %============Start computations========================
  
  nx = length(x);
  dx = x(2) - x(1);
  nx1 = round((x1-x(1))/dx)+1;  % x index corresponding to x1
  nx2 = round((x2-x(1))/dx)+1;  % x index corresponding to x2
  
  erfcx = erfc(x); % exact solution
  
  % Leading (n=0) term tn and partial sum sn=tn in series
  
  an = 0.5;
  tn = 2/sqrt(pi) * exp(-x.^2).*an./x;
  sn = tn;
  
  subplot(2,2,1)
  semilogy(x,erfcx,'b-',x,sn,'r--')  % Plot visual fit of leading-order term.
  xlabel('x')
  title('Large-x approximation to erfc(x)')
  legend('erfc(x)','t_0(x)')
  errnm1 = zeros(P+1,nx);  % Residual error matrix
  rat = zeros(P+1,nx);    % Matrix of ratios of n'th term to error
                          % after n-1 terms
  rat(1,:) = tn./erfcx;
  errnm1(1,:) = erfcx - sn;  % Residual before term 1

  for n = 1:P
    
    % Iteratively calculate next term and partial sum in series
    
    tn = -tn*(n-0.5)./x.^2;
    sn = sn + tn;
    rat(n+1,:) = tn./errnm1(n,:);
    errnm1(n+1,:) = erfcx - sn;
  end

  subplot(2,2,2)
  plot(x,rat(1:5,:))
  ylim([0 3])
  legend('n=1','2','3','4','5','Location','SouthWest')
  hold on
  plot(x,ones(1,nx),'k--',x,2*ones(1,nx),'k:')
  text(0.5,2.5,'t_n increases error')
  text(2,1.25,'t_n reduces error')
  hold off
  xlabel('x')
  ylabel('t_n(x)/\epsilon_{n-1}(x)')
  title('Ratio of new term to error')
  
  subplot(2,2,3)
  semilogy(0:P,abs(errnm1(:,nx1)))
  xlabel('P')
  ylabel('|\epsilon_P(x)|')
  title(['Error vs. P for x = ' num2str(x1)])
  
  subplot(2,2,4)
  semilogy(0:P,abs(errnm1(:,nx2)))
  xlabel('P')
  ylabel('|\epsilon_P(x)|')
  title(['Error vs. P for x = ' num2str(x2)])