clear; format long;
ns=400;

load sg_c05.dat
u=sg_c05(:,1)'+0e-06*rand(1,ns);
p=zeros(1,ns);

it=1; dt=0.002; cc=0.5; dx=sqrt(1/cc);
it5=1;it2=1;

while ((it-1)*dt<500)

k1u=dt*p; 
k1p=dt*(cc*[0,diff(u,2),0]-sin(u));

a=u+k1u/2;
b=p+k1p/2;

k2u=dt*b; 
k2p=dt*(cc*[0,diff(a,2),0]-sin(a));

a=u+k2u/2;
b=p+k2p/2;

k3u=dt*b; 
k3p=dt*(cc*[0,diff(a,2),0]-sin(a));

a=u+k3u;
b=p+k3p;

k4u=dt*b; 
k4p=dt*(cc*[0,diff(a,2),0]-sin(a));

u1=u+(k1u+2*k2u+2*k3u+k4u)/6;
p1=p+(k1p+2*k2p+2*k3p+k4p)/6;

if (mod(it,50)==0)

  tim(it2)=(it-1)*dt;
  probe(it2)=u(ns/2+1);     
     for i=1:ns-1
       if ((u1(i+1)>pi) & (u1(i)<pi))
	 xc=(i+(pi-u1(i))/(u1(i+1)-u1(i)));
       end;
     end;

pos(it2)=xc;

  it2=it2+1;

% the commands below allow visualization on the fly (i.e., as the code runs)  
%  subplot(2,1,1)
%plot(linspace(1,ns,ns)-ns/2,u,'-o')
%drawnow
%subplot(2,1,2)
%plot(tim,pos);
%drawnow;



if (mod(it,1000)==0)
u_save(:,it5)=u;
it5=it5+1
end;

end;

u=u1; p=p1;
it=it+1;
end;