eps=input('  eps=  ');

    betaL=1; betaR=2;

 for NN=1:4

    N=16*(2^NN);
    D1=DM1(N);  D2=D1^2;

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

    for i=1:N-1
       s=x(i); p1=-s;  q1=-1; 
       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);
       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';

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

       xx(NN)=N;
       err1(NN)=error;

 end
    fprintf(1,    '%10.0f        %10.3e \n', [xx; err1]);