%code 2.2 for lecture 2; need u_0(x), u_L(t), u_R(t)

    dt=input('  dt=  ');
    T=input(  'T=  ');
    epp=input(  'epp=  ');

    x1=1;  x2=-1;

for N=3:10
   
    j=[1:1:N-1];  x=[cos(j*pi/N)]';  U=epp*(x+tan(0.5*x));

    D1=DM1(N); D2=D1*D1;
    A=D2(2:N,2:N);  B=D1(2:N,2:N);

     time=0.0;

    while time < T 

       U0=U;

       C=diag(U);

       guo=1+epp*time;
       UR=(epp/guo)*(x1+tan(0.5*x1/guo));
       UL=(epp/guo)*(x2+tan(0.5*x2/guo));

       b=RK_RHSb(N,UL,UR,epp,D1,D2,U);

       K1=epp*A*U-C*B*U+b';

       U=U0+0.5*dt*K1;
       C=diag(U);
       
       guo=1+epp*(time+0.5*dt);
       UR=(epp/guo)*(x1+tan(0.5*x1/guo));
       UL=(epp/guo)*(x2+tan(0.5*x2/guo));
      
       b=RK_RHSb(N,UL,UR,epp,D1,D2,U);
 
       K2=epp*A*U-C*B*U+b';

       U=U0+0.5*dt*K2;
       C=diag(U);

       b=RK_RHSb(N,UL,UR,epp,D1,D2,U);
 
       K3=epp*A*U-C*B*U+b';

       U=U0+dt*K3;
       C=diag(U);

       guo=1+epp*(time+dt);
       UR=(epp/guo)*(x1+tan(0.5*x1/guo));
       UL=(epp/guo)*(x2+tan(0.5*x2/guo));
      
       b=RK_RHSb(N,UL,UR,epp,D1,D2,U);
 
       K4=epp*A*U-C*B*U+b';

       U=U0+dt*(K1+2*K2+2*K3+K4)/6;

       time=time+dt;

    end

%exact solution

       guo=1+epp*T;
       exact=(epp/guo)*(x+tan(0.5*x/guo));
       error=abs(exact-U);

       xx(N-2)=N;
       err(N-2)=max(error);
       ss=max(error)

end
  
    fprintf(1,      '%9.0f          %10.2e  \n', [xx;  err]);