#! /usr/bin/octave -qf

dummy = [0 : (nsteps)];
S = [dummy/sw];
U = [ceil(wr*dummy/(Xmax*n*sw))];
U(1)=1;
sol= zeros (nsteps+1, 2);
th = zeros (2,1);
global Xj
Xj=[S(dummy+1)*wr./(U(dummy+1)*n)];
global i
i = 1;

function csa = f(th)
global i 
global n
global Xj
csa(1)=(sin(2*pi*th(1))-sin(2*pi*th(2))+sin((n*pi/2)-th(1)*2*pi+th(2)*2*pi)-pi*n*Xj(i)/2) ;
csa(2)=(sin(4*pi*th(1))-sin(4*pi*th(2))+sin((n*pi)-th(1)*4*pi+th(2)*4*pi)-pi*n*Xj(i)) ;
endfunction 

for k = [1:nsteps+1];
clear th ;
i = k ; 
[th, info] = fsolve( "f", [(0.2*n/2) (0.4*n/2)]) ;
sol(k,:)=th ; 
endfor;
vals3=[n/4+sol(:,2)-sol(:,1)] ;
output = [Xj' (U-1)' sol vals3] 
vdlist = [sol vals3]' ;
vclist = [(U-1)'] ;