function [ar] = analytic_1_d_flux   (i_p, ax, Exact)
% Copyright 2000 J.E. Akin. All rights reserved.
%  variable "ar" (analytic flux  ) for "Exact" case number
%  of i_p-th component value at locations ax= [xmin:a_inc:xmax]

  if ( nargin < 3 )
    error ('No solution given for Exact_Case number')
  end % if 1 argument

  %b UNDER_CONSTRUCTION look for %OK
  if ( i_p >= 1 )  % actual component

    if ( Exact == 1 ) % U_ext + (U_0-U_ext) * Cosh (m*(L-x) / COSH_ML
      L     = input ('Enter physical length, L ') ;
      U_ext = input ('Enter convection temperature, U_ext ') ;
      U_0   = input ('Enter essential bc at x = 0, U_0 ') ;
      m     = input ('Enter m = sqrt (h_e*P_e/(K_e*A_e)) ') ;
      cosh_ML = cosh (m*L) ;
      ar = -(U_0-U_ext)*m*sinh(m*(L-ax))/cosh_ML; %OK  

    elseif ( Exact == 9 ) % u" + U + x = 0, EBC, EBC
      ar = cos(ax)/sin(1) - 1. ; % analytic flux  %OK

    elseif ( Exact == 10 ) % u" + U + x = 0, EBC, NBC
      ar = cos(ax)/cos(1) - 1. ; % analytic flux  %OK

    elseif ( Exact == 11 ) % u" + X^N = 0,  U(0)=0=U(1),
      N = input ('Enter source exponent, x^N ')
      ar = (1.-(N+2)*ax.^(N+1))/((N+1)*(N+2)) ; % analytic flux %OK

    elseif ( Exact == 14 ) % u" + Q_step = 0,  U(0)=0=U(1)
      for j = 1:(10*nt + 1) % loop over points
        if ( ax(j) <= 0.5 )
          ar(j) = (3 - ax(j))/8 ; %OK
        else
          ar(j) = -1. /8 ; %OK
        end % if position
      end % for

    elseif ( Exact == 16 ) 
      % Y'' - 2XY'/(X^2+1) + 2Y/(X^2+1) = (X^2+1), Fausett p. 484
      % with Y(0)=1, Y'(1)+Y(1)=0, Y=(X^4 -3X^2 -X +6)/6
      ar = (4* ax.^4 - 6*ax -1.)/6. ; %OK

    elseif ( Exact == 28 ) % u dp/dx - d(1 * dp/dx)/dx = 3u*x^2
      U = 1 ; % 60 ;
      const = exp(U) - 1;
      for j = 1:(10*nt + 1) % loop over points
        %b ar(j) =  ax(j)^3 + 3*ax(j)^2/U + 6*ax(j)/U/U   ...
        %b       - (1 + 3/U + 6/U/U)*(exp(U*ax(j))-1)/const ;
        ar(j) = 0. ; % BAD
      end % for

    elseif ( Exact == 29 ) % u dp/dx - d(1 * dp/dx)/dx = 3u*x^2 +
                   % u*TWO_PI*(COS(TWO_PI*X)) + TWO_PI*SIN(TWO_PI*X))
      U = 60 ; % 1 ; % 60 ;
      const = exp(U) - 1;
      for j = 1:(10*nt + 1) % loop over points
        %b ar(j) =  ax(j)^3 + 3*ax(j)^2/U + 6*ax(j)/U/U   ...
        %b       - (1 + 3/U + 6/U/U)*(exp(U*ax(j))-1)/const ...
        %b       + sin (2*pi*ax(j)) ;
        ar(j) = 0. ; % BAD
      end % for
    else
      error ('No solution given for Exact_Case number')
    end % if Exact

  %b else % i_p = 0, get root mean sq
  end % if i_p >= 1

  if ( i_p == 0 ) % get root mean sq

    if ( Exact == 9 ) % u" + U + x = 0, EBC, EBC
      ar = abs(cos(ax)/sin(1) - 1.) ; % analytic flux  %OK

    elseif ( Exact == 10 ) % u" + U + x = 0, EBC, NBC
      ar = abs(cos(ax)/cos(1) - 1.) ; % analytic flux  %OK

    elseif ( Exact == 11 ) % u" + X^N = 0,  U(0)=0=U(1),
      N = input ('Enter source exponent, x^N ')
      ar = abs(1.-(N+2)*ax.^(N+1))/((N+1)*(N+2)) ; % flux %OK

    elseif ( Exact == 16 ) 
      % Y'' - 2XY'/(X^2+1) + 2Y/(X^2+1) = (X^2+1), Fausett p. 484
      % with Y(0)=1, Y'(1)+Y(1)=0, Y=(X^4 -3X^2 -X +6)/6
      ar = (4* ax.^4 - 6*ax -1.)/6. ; %OK
      ar = abs (ar) ;

    else
      error ('No solution given for Exact_Case number')
    end % if Exact

  end % if get RMS value

% end function analytic_1_d_flux