function [k,m] = myTISE_C1_L2(x_e,pt,wt)   % SHO Schrodinger element
%------------------------------------------------------------------- 
%  Purpose: 
%     TISE stiffness and mass matrices for Hermitian line element 
%     nodal dof {v_1 theta_1 v_2 theta_2} , theta = v,x
% 
%  Variable Description: 
%     k - element stiffness matrix (size of 4x4)    
%     j - element mass matrix (size of 4x4) 
%     p - by exact numerical integration
%     h - current element length 
%     n - element number
%     pt,wt - Gaussian quadrature data, on -1 to +1
%     x - x at quadrature point    
%------------------------------------------------------------------- 
 
 h = x_e(2) - x_e(1);     % physical length
 
% exact consistent mass matrix: Integral, over h, Y' * Y dx
% ' is transpose (complex conjugant)
 
  m = h * [156      22*h     54     -13*h  ;... 
           22*h     4*h^2    13*h   -3*h^2 ;... 
           54       13*h     156    -22*h  ;... 
          -13*h    -3*h^2   -22*h    4*h^2] / 420; 
 
  k   = zeros (4, 4); % initialize k
  p   = zeros (4, 4); % initialize p
  Jac = h / 2;     % Jacobian for quadrature rule
  nqp = size(pt,2);   % rule number (here 5 will be exact)

  for j = 1:nqp       % loop over quadrature points
    r = pt(j);        % non-dimensional position
    x = ((1-r)*x_e(1) + (1+r)*x_e(2))/2; % physical x
    pot = x^2;        % SHO
    %b [pot] = potential (x,lib);   % User defined potential function

% form C1 cubic interpolation functions, Y = N_C1*Y_nodes
    N_C1(1) = (2 - 3*r + r^3)/4;          
    N_C1(2) = (1 - r - r^2 + r^3)*h/8;
    N_C1(3) = (2 + 3*r - r^3)/4;    
    N_C1(4) = (-1 - r + r^2 + r^3)*h/8;

% potential matrix: p = Integral, over h, Pot(x) * Y' * Y dx
    p = p + pot * N_C1' * N_C1 * Jac * wt(j) ;

% form slope of C1 cubic interpolation functions, Y,x = DN_C1*Y_nodes
    DN_C1(1) = (-3 + 3*r*r) * 0.5 / h;
    DN_C1(2) = (-1 - 2*r + 3*r*r) * 0.25 ;
    DN_C1(3) = ( 3 - 3*r*r) * 0.5d0 / h;
    DN_C1(4) = (-1 + 2*r + 3*r*r) * 0.25;
 
% diffusion matrix: Integral, over h,  Y',x * Y,x dx 
    k = k + DN_C1' * DN_C1 * Jac * wt(j) ;

  end % for j numerical integration loop

%    add potential to diffusion term
  k = k + p;    % final stiffness

%b  Second derivatives for later Ham, <Y|Ham|Y>, etc
%b D2N_C1(1) = ( 6*r) / (h*h)
%b D2N_C1(2) = (-2 + 6*r) * 0.5d0 / h
%b D2N_C1(3) = (-6*r) / (h*h)
%b D2N_C1(4) = ( 2 + 6*r) * 0.5d0 / h

% end function myTISE_C1_L2