eps=input('  eps=  ');

    betaL=-4*exp(-5)+sin(1)+2*cos(1);  
    betaR=6*exp(5)+sin(1)+2*cos(1);
    bL=1; BL=-1; bR=1; BR=1;

for N=5:30

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

    ta=BL*D1(N+1,1);  tb=bL+BL*D1(N+1,N+1);
    tc=bR+BR*D1(1,1);  td=BR*D1(1,N+1);
    te=ta*td-tc*tb;

    tbR=(td*betaL-tb*betaR)/te;
    tbL=(ta*betaR-tc*betaL)/te;

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

    for i=1:N-1
       s=x(i); p1=s;  q1=-1;
       u_true(i)=exp(5*s)+sin(s^2);
       f1=(24+5*s)*exp(5*s)+(2+2*s^2)*cos(s^2)-(4*s^2+1)*sin(s^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
          ss3=(td*BL*D1(N+1,j+1)-tb*BR*D1(1,j+1))/te;
          ss4=(ta*BR*D1(1,j+1)-tc*BL*D1(N+1,j+1))/te;
          ss5=ss1*ss3+ss2*ss4;
            if i==j
              A(i,j)=eps*D2(i+1,j+1)+p1*D1(i+1,j+1)+q1-ss5;
            else
              A(i,j)=eps*D2(i+1,j+1)+p1*D1(i+1,j+1)-ss5;
            end
          end
       b(i)=f1-ss1*tbR-ss2*tbL;
    end

    u=A\b';

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

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

    end

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