function [PHI, DPHI, DDPHI] = MLS1DShape(m, nnodes, xi, npoints, x, dm, wtype, para)
wi  = zeros (1, nnodes);  % Weight funciton
dwi  = zeros (1, nnodes);
ddwi = zeros (1, nnodes);

PHI  = zeros(npoints, nnodes);
DPHI  = zeros(npoints, nnodes);
DDPHI = zeros(npoints, nnodes);

for j = 1 : npoints

for i = 1 : nnodes
di = x(j) - xi(i);
      [wi(i), dwi(i), ddwi(i)] = Weight(wtype, para, di, dm(i));
end
 
  if (m == 3)
      p = [ones(1, nnodes); xi; xi.*xi];
      px  = [1; x(j); x(j)*x(j)];
      dpx  = [0; 1; 2*x(j)];
      ddpx = [0; 0; 2];
     
      B    = p .* [wi; wi; wi];
      DB  = p .* [dwi; dwi; dwi];
      DDB  = p .* [ddwi; ddwi; ddwi];
  else
      error('Invalid order of basis.');
  end
 
A  = zeros (m, m);
DA  = zeros (m, m);
DDA = zeros (m, m);
for i = 1 : nnodes
      pp = p(:,i) * p(:,i)';
     
      A  = A  + wi(i) * pp;
      DA  = DA  + dwi(i) * pp;
      DDA = DDA + ddwi(i) * pp;
  end
 
  AInv = inv(A);
     
  rx  = AInv * px;
  PHI(j,:) = rx' * B;  % shape function
   
  drx  = AInv * (dpx -DA * rx);
  DPHI(j,:) = drx' * B + rx' * DB;  % first order derivatives of shape function
 
  ddrx  = AInv * (ddpx - 2 * DA * drx - DDA * rx);
  DDPHI(j,:) = ddrx' * B + 2 * drx' * DB + rx' * DDB;    % second order derivatives of shape function
end