% Emacs mode header: -*-Octave-*- source. %
tic
1;
clear
N=10000;
global L=10;
global alpha=1;
global omega0=1000;
global d=2;
global omega=20;
global gamma=-0.0;
global B=1e-7;
rho=1;
global mu=rho*omega*omega;
function zdot = func(z,x0)
     global mu;
     global d;
     global gamma;
     global B;
     A=lin_res_A(x0);
     p2=z(1)*z(1)+z(2)*z(2);
     u=u_nlin(x0,p2);
     Lr = -real(mu/(A+i*B*u+i*gamma));
     Li = -imag(mu/(A+i*B*u+i*gamma));
     zdot(1)=z(3);
     zdot(2)=z(4);
     zdot(3)=Lr*z(1)-Li*z(2);
     zdot(4)=+Lr*z(2)+Li*z(1);
endfunction

function A = lin_res_A(x0)
     global alpha;
     global omega0;
     global L;
     global omega;
     global d;
     A=alpha*(omega0*exp(-(L-x0)/d)-omega);
endfunction 

function u = u_nlin(x0,p2)
     global B;
     global gamma;
     if (B==0)
      u=0;
     else
        a0=-p2/B^2;
        a1=(lin_res_A(x0)^2+gamma^2)/(B^2);
        a2=2*gamma/B;
     	s=a1/3-a2*a2/9;
     	r=-a0/2+a1*a2/6-a2^3/27;
        u = (r+sqrt(s^3+r^2))^(1/3)+...
        sign(r-sqrt(s^3+r^2))*(abs(r-sqrt(s^3+r^2)))^(1/3);
     endif
endfunction

npoint=9;
p0=[1e-9,5e-9,1e-8,5e-8,1e-7,5e-7,1e-6,5e-6,1e-5];

for l=1:npoint

k0=omega*sqrt(-rho/(lin_res_A(0)));
q0r=real(k0); %
q0i=imag(k0); %

pabs0=p0(l); %
z0=[0,0,pabs0*q0r,pabs0*q0i];
x=linspace(0,L,N);
zres=lsode("func",z0,x);

p=zres(:,1)+i*zres(:,2);
q=zres(:,3)+i*zres(:,4);

for j=1:N
 p2=abs(p(j)*p(j));
 x0=x(j);
 h(j)=p(j)/(lin_res_A(x0)+i*B*u_nlin(x0,p2)+i*gamma);
 u(j)=u_nlin(x0,p2);
endfor

xr=L-d*log(omega0/omega);
n=(N-1)*(xr/L)+1;
%p20(l)=abs(p(N)*p(N));
%p2x(l)=abs(p(n)*p(n));
%h2x(l)=abs(h(n)*h(n));

rh=real(h);
ih=imag(h);

% gset nologscale;
% data = [x',rh];
% gplot data w l

%filename=["data",num2str(l),".dat"];
%myfile=fopen(filename,"w");
%fprintf(myfile,"%10.5e %10.5e \n",rh,ih);
%fclose(myfile);
endfor

toc
%data=[p20,h2x];

%gset logscale xy
%gplot data w line
%%% p20
%%% h2x
