function T82
N=200;
eps=1/10;
a=2 ;
x=linspace(0,2*pi,N)';
r0=(x(1)+x(end))/2;
x0=r0;
h=x(2)-x(1);
df=@(r)1./(power(eps,2) + power(-r0 + r,2));
mean_df=@(r)df(r)./sqrt(df(r).^2+1);
d1u=df(x);
ss=@(x)((-2*x + 2*x0)./(sqrt(1 + power(power(eps,2) + power(x - x0,2),-2)).*...
     (1+power(eps,4)+power(x,4) + 2*power(eps,2).*power(x - x0,2) - 4*power(x,3)*x0 +... 
       6*power(x,2).*power(x0,2) - 4*x.*power(x0,3) + power(x0,4))));
An=-(atan((r0 - x)/eps)/eps);
S1=-ss(x);
D_BC=[mean_df(x(1:2)-2*h),mean_df(x(end-1:end)+2*h)];
N_BC=[df(x(1)),df(x(end))];

u=sin(x);
t=0;
dt=2*h;
t_end=2000;
L=eye(N,N)*(1/dt)-a*D2(x);
%--------------------------------------------
while t<t_end
     f_old=L\(S(u));
     u=f_old;
     t=t+dt;
     plot(x,u,'b-.',x,An,'r')
     title(['t=',num2str(t)])
    
     drawnow
end


%--------------------------------------------
function s=S(u)
         du=D3(x,u,N_BC,'N');
         p=du./sqrt(1+d1u.^2);
        %-----2. 外部使用WENO5---------
         d2u_p=WENO3(x,p,1,'no');
         d2u_m=WENO3(x,p,-1,'no');
         J=(d2u_p+d2u_m)/2;
         s=u/dt+J+S1-a*d2u(x,u);
    end
%--------------------------------------------
function y=d2u(x,u)
        % The goal of this function is calculate the second order
        % derivative  of  u, with a periodic boundary condition
          h=x(2)-x(1);
        um1=circshift(u,[1,0]);
        up1=circshift(u,[-1,0]);
        up1(end)=2*h*N_BC(2)*0 + u(end-1);
        um1(1)  =-2*h*N_BC(1)*0 + u(2);
        y=(um1-2*u+up1)/(h^2);
    end

    function y=D2(x)
        % out put: Matrix D2, with periodic boundary condition
         h=x(2)-x(1);
         n=length(x);
         c1=ones(n-2,1);
         c2=ones(n-1,1);
         c3=-2*ones(n,1);
         M=sparse(diag(c2,-1)+diag(c3,0)+diag(c2,1));
         M=M+sparse([1,n],[2,n-1],[1,1],n,n);
         y=M/(h^2);
    end
end