function neurls(N,M)
% N : nb of neurons
% M : nb of interpolation points
% mat: the rectangular matrix corresponding to the least square problem
% val: the values at interpolation points
% weight : the weights of the neural network
  xinterp=zeros(M,1);
  val=zeros(M,1);
  for j=1:M
    xinterp(j)=j/(M+1);
    val(j)=1+tanh(100*(xinterp(j)-0.5));
  end
  mat=zeros(M,N);
  for i=1:N
    for j=1:M
      mat(j,i) = basis(N,i,???);
    end	
  end
  weight=(mat'*mat)\(mat'*val)
% plot results
  xx=zeros(100,1);
  yy=zeros(100,1);
  for k=1:100
    xx(k)=k/101;
    yy(k)=sumbasis(N,weight,xx(k));
  end
  axis([0 1 -1 1]);
  plot(xinterp,val,'o',xx,yy);
  refresh();
end
function sumbasis=sumbasis(N,weight,x)
  sumbasis=0;
  for i=1:N
    sumbasis=sumbasis+weight(i)*basis(N,i,x);
  end
end
function basis=basis(N,i,x)
  basis=max(0,x-i/(N+1));
end

