%--------------------------                            
% program winkler.m
%--------------------------                            
 clear

%---------------- 
% INPUT DATA
%---------------- 
  load HermiteMESH  -ASCII
  load HermiteQUAD   -ASCII

% ---------------- 
% Define HermiteMESH Data
% ---------------- 
  NUMNP = HermiteMESH(1,1); 
  NUMEL = NUMNP-1;
  NUMEQ = 2*NUMNP;
  IB    = 4;
  for i=1:NUMNP;
     XORD(i)   = HermiteMESH(i+1,1);
     NPcode(i) = HermiteMESH(i+1,2);
  end

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

% --------------------------- 
% General initialization
% --------------------------- 
  for I=1:NUMEQ
     LHS(I)=0.0;
     RHS(I)=0.0;
     for J=1:4
        SK(I,J)=0.0; 
     end
  end
  for I=1:NUMNP
     NPBC(I)=0.0;
  end

  HermiteINITIAL

   NPBC

% -------------------------------- 
% FORMATION OF STIFFNESS MATRIX   
% and RIGHT HAND SIDE
% -------------------------------- 

  for I=1:NUMEL
% ----------------------------
% Initialize Element Variables
% ----------------------------
    for J=1:4
      for K=1:4
         S1(J,K)=0.0; 
         S2(J,K)=0.0; 
       end
       Qe(J)=0.0; 
    end

% ----------------------------
% Calculate element coordinate
% information.
% ----------------------------
     Xa = XORD(I);
     Xb = XORD(I+1);
     RL=Xb-Xa; 
     DxDu=RL/2.0; 
     DuDx=2.0/RL;   

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

%    ----------------------------
%    INCLUDE COEF.m 
%    Defines: EIx, Kx, and Qx
%    ----------------------------
     HermiteCOEF

%    ------------------------------
%    Calculate  shape functions wrt 
%    x at current Gauss point
%    ------------------------------
     No(1) = HermiteSF(0,1,u);
     No(2) = HermiteSF(0,2,u)*(DxDu);
     No(3) = HermiteSF(0,3,u);
     No(4) = HermiteSF(0,4,u)*(DxDu);

     Nx(1) = HermiteSF(1,1,u)*(DuDx);
     Nx(2) = HermiteSF(1,2,u)*(DuDx)*(DxDu);
     Nx(3) = HermiteSF(1,3,u)*(DuDx);
     Nx(4) = HermiteSF(1,4,u)*(DuDx)*(DxDu);


% ----------------------------
% Element stiffness matrices
% ----------------------------
    Wt = GWTS(J);
    for K=1:4
      for L=1:4
        S1(K,L)=S1(K,L) + ...
          Wt*Nx(K)*EIx*Nx(L)*DxDu;

        S2(K,L)=S2(K,L) + ...
          Wt*No(K)*Kx*No(L)*DxDu;
      end   
      Qe(K)=Qe(K) + Wt*No(K)*Qx*DxDu; 
   end   
  end   
% --------- Quadrature now complete 


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


% -------------------------
% BOUNDARY CONDITIONS
% -------------------------
  B = 1.0E+06;
  for I=1:NUMNP
   if NPBC(I) == 1 
     I1 = 2*I-1; 
     SK(I1,1)=SK(I1,1)*B;
     RHS(I1)=LHS(I1)*SK(I1,1);
   end
  end


% -------------------------
% CALL EQUATION SOLVER
% -------------------------
  LHS = sGAUSS(SK,RHS,NUMEQ,IB);

% --------------------------- 
% Solution is now complete.
% Save values.
% --------------------------- 

  save LHS LHS -ASCII

      figure(2)
      hold on;
      for I=1:NUMEL
          L_I=XORD(I+1)-XORD(I);
          h_I= L_I/10;
          x_I= XORD(I):h_I:XORD(I+1);
          u_I= (2*x_I-XORD(I+1)-XORD(I))/L_I;
          y_I =LHS(2*I-1)* HermiteSF(0,1,u_I)+LHS(2*I)* HermiteSF(0,2,u_I)*L_I/2 ...
              +LHS(2*I+1)* HermiteSF(0,3,u_I)+LHS(2*I+2)* HermiteSF(0,4,u_I)*L_I/2;
          app=plot(x_I,y_I,'--r', 'linewidth', 2);
          grid on;
      end
      title(' Results from HermiteOde2.m')
      xlabel(' Independent variable')
      ylabel(' Dependent variable')
     
      x1 = min(XORD):0.01:0;
      x2 = 0:0.01:max(XORD);
      %exact=plot(x1,x1, '-k', 'linewidth', 2);
      %plot(x2,x2, '-k', 'linewidth', 2);
      a=20
      exact=plot(x1,2*x1/(a+1)+(1-a)/(a+1), '-k', 'linewidth', 2);
      plot(x2,2*a*x2/(a+1)+(1-a)/(a+1), '-k','linewidth', 2);
      legend([app exact], 'Approximate solution', 'Exact solution');
      
      hold off
      


