%===============
%  PROGRAM ode2
%===============

      clear
%     ---------------- 
%     INPUT DATA
%     ---------------- 
      load MESH   -ASCII
      load QUAD   -ASCII

%     ---------------- 
%     Define MESH Data
%     ---------------- 
      NUMNP = MESH(1,1); 
      NUMEL = NUMNP-1;
      for i=1:NUMNP;
        XORD(i) = MESH(i+1,1);
        NPBC(i) = MESH(i+1,2);
        Y(i)    = MESH(i+1,3);
        Q(i)    = 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

%     --------------------------- 
%     Initialize stiffness matrix
%     --------------------------- 
      for I=1:NUMNP
        for J=1:3
          SK(I,J)=0.0; 
        end
      end

%     -------------------------------- 
%     FORMATION OF STIFFNESS MATRIX   
%      and 
%     RIGHT HAND SIDE
%     -------------------------------- 
      IB=3;
      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;   

        for J=1:2
          for K=1:2
            S1(J,K)=0.0; 
            S2(J,K)=0.0; 
            S3(J,K)=0.0; 
          end
          Qe(J)=0.0; 
        end

%       ------------------------- 
%       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

%         ------------------------------
%         Gradient of shape functions at 
%         current Gauss point
%         ------------------------------
          DNDX(1) = SF(1,1,u)*DuDx;
          DNDX(2) = SF(1,2,u)*DuDx;
 
%         ----------------------------
%         Element stiffness matrices
%         ----------------------------
          Wt = GWTS(J);
          for K=1:2
            for L=1:2
              S1(K,L)=S1(K,L) +Wt*DNDX(K)*AX*DNDX(L)*DxDu;
              S2(K,L)=S2(K,L) +Wt*SF(0,K,u)*BX*DNDX(L)*DxDu;
              S3(K,L)=S3(K,L) +Wt*SF(0,K,u)*CX*SF(0,L,u)*DxDu;
            end   
            Qe(K)=Qe(K) + Wt*SF(0,K,u)*DX*DxDu; 
          end   
        end   
%       --------- Quadrature now complete 

%       --------------------------
%       Assemble element matrices 
%       into global matrix 
%       -------------------------- 
        K1=I-1;
        L0=2;
        for K=1:2
          L0=L0-1; 
          K1=K1+1; 
          L1=L0; 
          for L=1:2
            L1=L1+1; 
            SK(K1,L1)=SK(K1,L1) +S1(K,L)-S2(K,L)-S3(K,L); 
          end
          Q(K1) =  Q(K1) + Qe(K);
        end
      end
%     ---- Global Matrices are assembled

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

%     -------------------------
%     CALL EQUATION SOLVER
%     -------------------------
      Y =  nGAUSS(SK,Q,NUMNP,IB);

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