eps=input('eps=  ');
    M=input('M=  ');

    betaL=1; betaR=2; pi1=pi/2;

for NN=1:5

    N=2^(NN+4)

    D1=DDM1(N);  D2=D1^2;

    j=[1:1:N-1];  y=[cos(pi*j/N)]';  

    for j=1:N-1
       gm=y(j);  yp(j)=1;
          for mm=1:M
             yp(j)=yp(j)*pi1*cos(pi1*gm);
             gm=sin(pi1*gm);
          end
        x(j)=gm;

       hm=x(j);  hd(j)=0;
          for mm=1:M
             hm=asin(hm)/pi1;
             hd(j)=pi1*tan(pi1*hm)+pi1*cos(pi1*hm)*hd(j);
          end
    end
      
    for i=1:N-1
       P1=-x(i)*yp(i)+eps*hd(i);
       Q1=-yp(i)^2;
       
       s=x(i); 
       s1=(s+1)/eps;  s2=(s-1)/eps;
       u_true(i)=exp(-s1)+2*exp(s2);

       f1=(s1-1)*exp(-s1)-2*(s2+1)*exp(s2);
       F1=f1*yp(i)^2;

       ss1=eps*D2(i+1,1)+P1*D1(i+1,1);
       ss2=eps*D2(i+1,N+1)+P1*D1(i+1,N+1);

          for j=1:N-1
            if i==j
              A(i,j)=eps*D2(i+1,j+1)+P1*D1(i+1,j+1)+Q1;
            else
              A(i,j)=eps*D2(i+1,j+1)+P1*D1(i+1,j+1);
            end
          end
       b(i)=F1-ss1*betaR-ss2*betaL;
    end

    u=A\b';

    for i=1:N-1
       error(i)=abs(u(i)-u_true(i));
    end

    xx(NN)=N;
    smx=max(error)
    err(NN)=smx;

end

    fprintf(1,      '%16.0f          %13.3e  \n', [xx;  err]);