eps=input('  eps=  ');

    betaL=exp(-5)+sin(1);  betaR=exp(5)+sin(1);

for N=5:30

    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;
       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
            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(N)=N;
    err(N)=max(error);

    disp(error)

    end

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