eps=input('  eps=  ');

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

for NN=1:6

    N=2^(NN+3)

    h=2/N; 

    j=[1:1:N-1];  x=-1+j*h;

    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);

          for j=1:N-1
            if i==j
              A(i,j)=-2*eps+q1*h^2;
            elseif i==j+1
              A(i,j)=eps-0.5*h*p1;
            elseif j==i+1
              A(i,j)=eps+0.5*h*p1;
            else
              A(i,j)=0;
            end
          end
      b(i)=f1*h^2;
    end

       ss1=eps+0.5*h*x(N-1);
       b(N-1)=b(N-1)-ss1*betaR;

       ss2=eps-0.5*h*x(1);
       b(1)=b(1)-ss2*betaL;

    u=A\b';

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

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

    end

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