%
%  f_loop.m
%     - modified to be a single file (no more fourthN.m)
%  			converged result;  K = -0.84193000177728
%	  - to start, set K=0 -> updates to K are now kept!
%
% - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 

%  timestep, steps/cycle & cycles & parameter
rkDt   = 0.005;  Nsteps = 100;  Ncycle = 36;
%rkDt   = 0.01;  Nsteps = 25;  Ncycle = 30;  

par = [1];  check = 1.0;

if(K==0);  K = -0.9;  end

while (abs(check)>1e-8)

%  initialize ode solve
	rkt = 0.0;  IVs = [0.0 K 0.0 0.0 1.0 0.0];  
	rkZ = IVs;  rkS = [rkt rkZ];

%  timestep loop
	for k = 1:Ncycle
  		for j = 1:Nsteps
			[rkZ,rkt] = fourth_rkN(rkZ,rkt,rkDt,par);
  		end
  
  		rkS = [rkS ; [rkt rkZ]];
	end

%  newton checks & update
	check  = rkS(end,2)+1 + rkS(end,3) + rkS(end,4);
	checkP = rkS(end,5)   + rkS(end,6) + rkS(end,7);
  
	K = K - check/checkP;
	disp(['K='num2str(K) ... 
			'  ;  error='num2str(check) ...
			'  ;  dk='num2str(-check/checkP)])
end

time  = [-flipud(rkS(2:end,1)) ; rkS(:,1)];
wave0 = [-flipud(rkS(2:end,2)) ; rkS(:,2)];

%  plot

plot(time,wave0,'k');  hold on
plot(time,wave0,'r.','markersize',7)

title('\bf hyper-diffusion front')
ylabel('\bf W(y)-axis');  xlabel('\bf y-axis')
pmax = max(wave0)*1.1;  pmin = -max(-wave0)*1.1;
axis([time(1) time(end) pmin pmax]);  hold off