clear
!clear
format long
%load fracna.dat

ns=200;
h1=(0.15);
cp=1/(2*h1^2);
h=0.001;
ii=1;

sp=linspace(1,ns,ns)*h1;

uin=sech(sp-(ns*h1)/2).*exp(i*0.2*(sp-ns*h1/2));
x=real(uin); y=imag(uin);
%+0.001*cos(2*pi*sp/3);
%x=exp(-(sp-(ns*h1)/2).^2/2);

it5=1;

%y=zeros(1,ns);

eps=0.0;
fl=1.0;
num=5;


gam=zeros(1,ns);
ep=1;
gam5=zeros(1,ns);
it=1;

while (ii<40000)


     d2y=diff(y,2);     ad1y=(y(2)+y(ns)-2*y(1)); ad3y=(y(1)+y(ns-1)-2*y(ns));
     d2x=diff(x,2);     ad1x=(x(2)+x(ns)-2*x(1)); ad3x=(x(1)+x(ns-1)-2*x(ns));
     p1=[ad1y,d2y,ad3y]; p2=[ad1x,d2x,ad3x];
     k1x=h*(-cp*p1-ep*(x.^2+y.^2).*y-gam.*y-gam5.*x);
     k1y=h*(cp*p2+ep*(x.^2+y.^2).*x+gam.*x-gam5.*y);
     a=x+k1x/2;
     b=y+k1y/2;

     d2y=diff(b,2);     ad1y=(b(2)+b(ns)-2*b(1)); ad3y=(b(1)+b(ns-1)-2*b(ns));
     d2x=diff(a,2);     ad1x=(a(2)+a(ns)-2*a(1)); ad3x=(a(1)+a(ns-1)-2*a(ns));
     p1=[ad1y,d2y,ad3y]; p2=[ad1x,d2x,ad3x];
     k2x=h*(-cp*p1-ep*(a.^2+b.^2).*b-gam.*b-gam5.*a);
     k2y=h*(cp*p2+ep*(a.^2+b.^2).*a+gam.*a-gam5.*b);
     a=x+k2x/2;
     b=y+k2y/2;

     d2y=diff(b,2);     ad1y=(b(2)+b(ns)-2*b(1)); ad3y=(b(1)+b(ns-1)-2*b(ns));
     d2x=diff(a,2);     ad1x=(a(2)+a(ns)-2*a(1)); ad3x=(a(1)+a(ns-1)-2*a(ns));
     p1=[ad1y,d2y,ad3y]; p2=[ad1x,d2x,ad3x];
     k3x=h*(-cp*p1-ep*(a.^2+b.^2).*b-gam.*b-gam5.*a);
     k3y=h*(cp*p2+ep*(a.^2+b.^2).*a+gam.*a-gam5.*b);
     a=x+k3x;
     b=y+k3y;

     d2y=diff(b,2);     ad1y=(b(2)+b(ns)-2*b(1)); ad3y=(b(1)+b(ns-1)-2*b(ns));
     d2x=diff(a,2);     ad1x=(a(2)+a(ns)-2*a(1)); ad3x=(a(1)+a(ns-1)-2*a(ns));
     p1=[ad1y,d2y,ad3y]; p2=[ad1x,d2x,ad3x];
     k4x=h*(-cp*p1-ep*(a.^2+b.^2).*b-gam.*b-gam5.*a);
     k4y=h*(cp*p2+ep*(a.^2+b.^2).*a+gam.*a-gam5.*b);

x1=x+(k1x+2*k2x+2*k3x+k4x)/6;
y1=y+(k1y+2*k2y+2*k3y+k4y)/6;
uu=x1.^2+y1.^2;

if (mod(ii,50)==0) 
tim(it)=ii*h;

pos2(it)=max(uu);
[s13,s14]=max(uu);
pos3(it)=sum(uu);
it=it+1;

xxxx=linspace(1,ns,ns);
subplot(2,1,1)
plot(xxxx,x1,'-',xxxx,y1,'--')
subplot(2,2,3)
plot(xxxx,uu)
drawnow
subplot(2,2,4)
plot(tim,pos3);
end;

if (mod(ii,200)==0) 
uu10(:,it5)=uu';
x10(:,it5)=x';
y10(:,it5)=y';
it5=it5+1;
end;
  
x=x1;
y=y1;

ii=ii+1;

end;