function [r] = Fcn2(x,p) r.x0=x(1,1); r.Kappa = x(2:p.n+1,1); p.y0=p.A*(sin(2*pi*p.t)); r.Theta = p.Theta0+cumtrapz(p.Grid, r.Kappa); r.Zeta = r.x0+1i*p.y0+cumtrapz(p.Grid, exp(1i*r.Theta)); r.hs = exp(1i*p.theta); %%%%%%\hat{s}%%%%%%%%%% r.hn=1i*r.hs; r.DtZeta=p.DtZeta; r.Dtx0=p.Dtx0; r.Zeta1=r.Zeta-2*p.dt*r.DtZeta; r.x01=r.x0-2*p.dt*r.Dtx0; r.DttZeta=(r.Zeta-2*p.Zeta+r.Zeta1)/p.dt^2; r.Dttx0=(r.x0-2*p.x0+r.x01)/p.dt^2; r.us=real(r.DtZeta.*conj(r.hs)); r.un=real(r.DtZeta.*conj(r.hn)); if (nnz([p.mub,p.muf,p.mut]) ~=0) NormDtZeta = abs(r.DtZeta); diff(1:p.n-1)=r.us(2:p.n,1)-r.us(1:p.n-1,1); a=diff/p.dx; b=r.us(2:p.n,1)'-p.Grid(2:p.n).*a; diff2(1:p.n-1)=r.un(2:p.n,1)-r.un(1:p.n-1,1); c=diff2/p.dx; d=r.un(2:p.n,1)'-p.Grid(2:p.n).*c; a2=2*(a.*b+c.*d)./(a.^2+c.^2); a3=(b.^2+d.^2)./(a.^2+c.^2); I0=log((a2+2*p.Grid(2:p.n)+2*sqrt(a3+a2.*p.Grid(2:p.n)+p.Grid(2:p.n).^2))./(a2+2*p.Grid(1:p.n-1)+2*sqrt(a3+a2.*p.Grid(1:p.n-1)+p.Grid(1:p.n-1).^2))); I1=sqrt(a3+a2.*p.Grid(2:p.n)+p.Grid(2:p.n).^2)-sqrt(a3+a2.*p.Grid(1:p.n-1)+p.Grid(1:p.n-1).^2)-0.5*a2.*I0; ind=find(r.us>0); Mufb=p.mub*ones(p.n,1); Mufb(ind)=p.muf; intemuf=Mufb(1:p.n-1)'./sqrt(a.^2+c.^2).*(a.*I1+b.*I0); r.T =cumtrapz(p.Grid, real(r.DttZeta.*conj(r.hs))+p.damp*r.us); r.T(2:p.n,1)=r.T(2:p.n,1)+cumsum(intemuf'); left=find(r.us(1:p.n-1,1).*r.us(2:p.n,1)<0)+1; for i=1:length(left) index=left(i); if index==2 v0=r.us(index-1,1); v1=r.us(index,1); v2=r.us(index+1,1); n0=p.Grid(index-1); n1=p.Grid(index); n2=p.Grid(index+1); po=[v0-2*v1+v2,-(n1+n2)*v0+2*(n0+n2)*v1-(n0+n1)*v2,n1*n2*v0-2*n0*n2*v1+n0*n1*v2]; zp=roots(po); us0=zp(find(zp>n0-1e-8&zpn1-1e-8&zp