
    function d=RK_RHSb(N,UL,UR,epp,D1,D2,U)

    for j=1:N-1
      ss=epp*(D2(j+1,1)*UR+D2(j+1,N+1)*UL);
      d(j)=ss-U(j)*(D1(j+1,1)*UR+D1(j+1,N+1)*UL);
    end

