%===============
%  PROGRAM Ode2
%===============
      clear
%     ---------------- 
%     INPUT DATA
%     ---------------- 
      load MESH   -ASCII
      load QUAD   -ASCII

%     ---------------- 
%     Define MESH Data
%     ---------------- 
      NUMNP = MESH(1,1); 
      NUMEL = NUMNP-1;
%     ----------------------- 
%     Declare sparse matrices
%     ----------------------- 
      SK=sparse(NUMNP,NUMNP);

      for i=1:NUMNP;
        XORD(i) = MESH(i+1,1);
        NPBC(i) = MESH(i+1,2);
        Y(i,1)    = MESH(i+1,3);
        Q(i,1)    = MESH(i+1,4);
      end

%     ---------------- 
%     Define QUAD Data
%     ---------------- 
      NQPTS = QUAD(1,1);
      for i=1:NQPTS
        GPTS(i)=QUAD(i+1,1);
        GWTS(i)=QUAD(i+1,2);
      end

%     -------------------------------- 
%     FORMATION OF STIFFNESS MATRIX   
%      and 
%     RIGHT HAND SIDE
%     -------------------------------- 
      for I=1:NUMEL
%       ----------------------------
%       Initialize Element Variables
%       ----------------------------
        Xa = XORD(I);
        Xb = XORD(I+1);
        RL=Xb-Xa; 
        DxDu=RL/2.0; 
        DuDx=2.0/RL;   

        S1(1:2,1:2)=0.0;
        S2(1:2,1:2)=0.0;
        S3(1:2,1:2)=0.0;
        Qe(1:2,1) = 0.0;

%       ------------------------- 
%       Begin Gaussian Quadrature 
%       ------------------------- 
        for J=1:NQPTS
          u = GPTS(J);
%         --------------------
%         Global coordinate of 
%         current Gauss point
%         --------------------
          Xg = (Xa+Xb)/2  + (RL/2.0)*u;

%         ----------------------------
%         INCLUDE COEF.m 
%         Defines: AX, BX, CX, and DX
%         ----------------------------
          COEF

%         ------------------------------
%         Shape functions and gradients
%         at current Gauss point and 
%         ------------------------------
          SFu(1)  = SF(0,1,u);
          SFu(2)  = SF(0,2,u);
          DNDX(1) = SF(1,1,u)*DuDx;
          DNDX(2) = SF(1,2,u)*DuDx;
 
%         ----------------------------
%         Element stiffness matrices
%         ----------------------------
          Wt = GWTS(J);
          S1 = S1+Wt*DNDX'*AX*DNDX*DxDu;
          S2 = S2+Wt*SFu'*BX*DNDX*DxDu;
          S3 = S3+Wt*SFu'*CX*SFu*DxDu;
          Qe = Qe+Wt*SFu'*DX*DxDu; 

        end   
%       --------- Quadrature now complete 

%       --------------------------
%       Assemble element matrices 
%       into global matrix 
%       -------------------------- 
        NP=[I,I+1];
        SK(NP,NP)=SK(NP,NP)+S1-S2-S3;
        Q(NP,1) =  Q(NP,1) + Qe;

      end
%     ---- Global Matrices are assembled

%     -------------------------
%     BOUNDARY CONDITIONS
%     -------------------------
      B = sum(sum(abs(SK)))*1.0E+04;
      for I=1:NUMNP
        if NPBC(I) == 1 
          SK(I,I)=B;
          Q(I)=Y(I)*B;
        end
      end

%     -------------------------
%     SOLVE EQUATIONS
%     -------------------------
      Y =  SK\Q;

%     ----------------------- 
%     OUTPUT DATA 
%     ----------------------- 
      save Y Y -ASCII;
      plot(XORD,Y)
      grid
      title(' Results from ode.m')
      xlabel(' Independent variable')
      ylabel(' Dependent variable')

