%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%This code finds type two explicit Multi-Step Runge-Kutta 
% Methods which preserves positivity
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%===============================================================
%Variable meanings:
% s         - # of stages
% p         - order of accuracy
% k         - # of steps
% starttype - type of initial guess (random, smart, load)

% Decision variables:
% A,d,b,theta - coefficients of the TSRK method
% r     - the radius of absolute monotonicity (multiplied by -1)
% x     - the decision variables stored in the following order:
% x=[A d' b' theta r] 
% A is stored row-by-row????
% minreff   - r/s
%===============================================================




% Fresh start:
%restart=0;  

rand('twister', sum(100*clock)); %New random seed every time
if restart==0 
 clear A D B Bhat Ahat theta
  restart=0;

  %Set optimization parameters:
%  opts=optimset('MaxFunEvals',1000000,'TolCon',1.e-15,'TolFun',1.e-15,'TolX',1.e-15,'GradObj','on','MaxIter',10000,'Diagnostics','on','Display','iter'...
%,'UseParallel','never','Algorithm','active-set');

opts=optimset('MaxFunEvals',100000,'TolCon',1.e-14,'TolFun',1.e-14,'TolX',1.e-14,'GradObj','on','MaxIter',10000,'Diagnostics','off','Display','off'...
,'UseParallel','never','Algorithm','sqp');
end

% Restart from known value
if restart~=0
opts=optimset('MaxIter',1000,'MaxFunEvals',10000,...
'TolCon',1.e-14,'TolFun',1.e-14,'TolX',1.e-14,...
'GradObj','on','Diagnostics','off','Display','off'...
,'Algorithm','sqp','UseParallel','never');
end


%==============================================
 %Editable options:
%  s=5; p=6; k=4;
   % solveorderconditions=1;
%  minreff = .1501; % Keep looking until method with at least this value is found
%==============================================

%==============================================
%Now set up:
%The number of unknowns -                     n
%The linear constraints -                     Aeq*x = beq
%The upper and lower bounds on the unknowns - ub, lb

n=.5*(s^2-s) + 2*s*k -s +1;
Aeq=[]; beq=[];
bd=20; % need to examine how large to pick bd
lb=zeros(1,n); lb(end)=-s*1.01;
ub=bd+zeros(1,n); ub(end)=0; 
%==============================================
count=0
info=-2;
while (info==-2 || (r/s)<minreff || info==0)

%while ((r/s)<minreff )
  %==============================================
  if restart==0
      if count==30
          ('exceed count')
          break
      end
  end
  
  if restart==1
    x=X;  
       if count==50
        ('exceed count')
        x=X1
        break
      end
  elseif restart==2
      if count==0
          X1=X;
      elseif count ~=0
          X=X1;
      end
      x=X+.2*rand(1,length(X));
      if count==50
        ('exceed count')
        x=X1
        break
      end
  
  else 
     x(1:n)=(rand(1,n));
	 x(n)=-0.01;
  
  end

  %==============================================
  %The optimization call:
  [X,FVAL,info]=fmincon(@tsrk_am_obj,x,[],[],Aeq,beq,lb,ub,@(x) nlc_tsrk(x,s,k,p,n),opts);
  r=-FVAL;
  
  count=count+1
  [r/s,info]
end %while loop
%==============================================
[A ,Ahat, B, Bhat, D, theta]  = unpackexplicitT2(X,s,k);

Aa= [zeros(k-1,s+k-1);Ahat,A];
Bb= [Bhat, B];
Dd= [eye(k);D];