%-----------------------------------------------------------------------
% Orthonormal functions. Numerical approach based on LU decomposition.
% Instructor: Nam Sun Wang
%-----------------------------------------------------------------------

% Start fresh ----------------------------------------------------------
      clear all

      global i j

% Program Header -------------------------------------------------------
      disp('Construct orthonormal functions e from the following given functions')
      disp('  f1(x)=J0(x) ')
      disp('  f2(x)=J0(2x)')
      disp('  f3(x)=J0(3x)')
      disp('  f4(x)=J0(4x)')
      disp('  f5(x)=J0(5x)')
      disp(' ')
      disp('The scalar product (f,g) is defined as')
      disp('          [1             ')
      disp('  (f,g) = [ f(x)*g(x)*dx ')
      disp('          [0             ')

% Find the scalar product between each pair of f
      n=5;
      for i=1:n
        for j=1:n
%         call a routine to find the integral
%         P(i,j)=quad8('orthogjf', 0., 1.); ... old
          P(i,j)=quadl('orthogjf', 0., 1.);
        end
      end

% Matlab's chol(P) returns the upper triangular matrix U such that U'*U=P
      U=chol(P);
      A=inv(U');

% Output results
      disp('Orthonormal basis functions e=A*f')
      disp('where A is...')
      disp(A)

% Plot orthonormal basis functions
      x=0. : 0.01 : 1.;
      i=0;
      y=A*orthogjf(x);
      plot(x,y)
      xlabel('x')
      ylabel('Orthogonal Basis Functions')
      axis([ 0. 1. -3. 3.])