function paramset(CI)
if nargin==0
%    CI=[0.05 0.06;0*0.05 0*0.06;0.2 0.23;0.05 0.055;0.3 0.35];%SEIRS
     CI=[0.2665681 0.2705561; 0.7157988 0.7586958; 0.0022428 0.0024325;1/5.3 1/5.3;1/5.61 1/5.61];%SIRCD
end

sCI=size(CI);
if sCI(2) ~= 2
    error('Confidence intervall must be a nx2 vector');
end

alpha=0.9;

%x0=[0.9;0;0.1;0];Iind=3;%SEIRS
x0=[10^6-1;0;1;0;1;0];Iind=6;%SEIRCD
p=(CI(:,2)+CI(:,1))./2;

Tfinal=300;%;273;
Tstep=Tfinal/100;
TotTime=(0:Tstep:Tfinal);

ObjVal=zeros(2,length(TotTime));
pextr=zeros(2,sCI(1),length(TotTime)-1);
Iextr=zeros(2,length(TotTime));
Iextr(:,1)=[x0(Iind) x0(Iind)];
for i=1:2
    l=zeros(length(x0),1);
    l(Iind)=(-1)^(i-1);
    ObjVal(i,1)=l'*x0;
    
    for j=2:length(TotTime)
        T=TotTime(j);
        
       Ttquot=2;
       time= (0:Tstep/Ttquot:T);
       [p,OV,state,costate] = grad_proj_loop(p,CI,l,@ffun,@fxfun,@fpfun,x0,time,alpha);
       Iextr(i,j)=state(Iind,end);
       ObjVal(i,j)=OV;
       pextr(i,:,j-1)=p;
       for k=1:j
           if (-1)^(i-1)*state(Iind,(k-1)*Ttquot+1)>(-1)^(i-1)*Iextr(i,k)
               Iextr(i,k)=state(Iind,(k-1)*Ttquot+1);
               pextr(i,:,k)=p;
           end
       end
        
       if 1==1
           figure(22)
           if i==1 && j==2
               clf;
               xlim([0 Tfinal]);
           end
           hold on
           if length(time)==2
               state=[state(:,1) state(:,end)];
           end
           plot(time,state(Iind,:)) 
       end
    end
end
 %%% for SEIRS
%figure(10);clf;
%plot(TotTime,Iextr(1,:))
%hold on
%plot(TotTime,Iextr(2,:))
%figure(11);clf;
%plot(TotTime,ObjVal(1,:))
%hold on
%plot(TotTime,ObjVal(2,:))
%figure(12);clf;
%for i=1:length(p)
%    subplot(length(p),2,2*i-1)
%    plot(TotTime(2:end),squeeze(pextr(1,i,:)))
%    if CI(i,1)~=CI(i,2)
%        ylim([CI(i,1) CI(i,2)])
%    end
%    subplot(length(p),2,2*i)
%    plot(TotTime(2:end),squeeze(pextr(2,i,:)))
%    if CI(i,1)~=CI(i,2)
%        ylim([CI(i,1) CI(i,2)])
%    end
%end


 %%%for SEIRCD
 meanstate=statesolv(@(t,x)ffun(t,x,[0.268574; 0.7368011;0.0023379;1/5.3;1/5.61]),TotTime,x0);%full
 figure(10)
 dataTime  =[111 113 114 115 116 117 120 121 124 127 129 131 134 136 137 140 143 151 153 155 157 160 162 168 173 177 178 182 184 186 196 197 201 211 213 217 219 223 225 228 231 234 238 241 243 246 248 251 253 255 258 260 262];
 dataCases =[49  86  86  86  103 112 122 127 143 151 158 159 167 197 203 208 218 226 231 235 236 233 248 253 258 281 291 328 344 351 394 398 390 413 412 408 409 406 411 410 415 427 460 472 485 495 495 506 510 519 543 579 607];
 dataDeaths=[29  59  60  62  66  70  80  83  86  95  101 106 108 122 129 136 141 149 155 157 158 157 171 176 174 186 193 209 215 226 263 264 270 303 305 307 309 304 310 310 314 319 339 346 358 363 367 373 377 380 394 396 406];
 if Iind==5
     scatter(dataTime(1:end),dataCases(1:end),'s')
     plot(TotTime,meanstate(5,:),'r:')
 elseif Iind==6
     scatter(dataTime(1:end),dataDeaths(1:end),'s')
     plot(TotTime,meanstate(6,:),'r:')
 end
 figure(923);clf;hold on
 meanstate=statesolv(@(t,x)ffun(t,x,[0.268574; CI(2,1); 0.0023379;1/5.3;1/5.61]),TotTime,x0);
 z=meanstate(5,:);
 meanstate=statesolv(@(t,x)ffun(t,x,[0.268574; CI(2,2); 0.0023379;1/5.3;1/5.61]),TotTime,x0);
 plot(TotTime,z-meanstate(5,:),'r')
end


%%SEIRS
%function dx=ffun(t,x,p)%nx1
%    dx(1,1)=-p(1).*x(1).*x(3) + p(2).*(x(2)+x(3)+x(4)) + p(3).*x(4);
%    dx(2,1)=p(1).*x(1).*x(3) - p(5).*x(2) - p(2).*x(2);
%    dx(3,1)=p(5).*x(2) - p(2).*x(3) - p(4).*x(3);
%    dx(4,1)=p(4).*x(3) - p(2).*x(4) - p(3).*x(4);
%end
%
%function fx=fxfun(t,x,p)%nXn
%    fx(1,1)=-p(1).*x(3);
%    fx(1,2)=p(2);
%    fx(1,3)=-p(1).*x(1)+p(2);
%    fx(1,4)=p(2)+p(3);
%    fx(2,1)=p(1).*x(3);
%    fx(2,2)=-p(2)-p(5);
%    fx(2,3)=p(1).*x(1);
%    fx(3,2)=p(5);
%    fx(3,3)=-p(2)-p(4);
%    fx(4,3)=p(4);
%    fx(4,4)=-p(2)-p(3);
%end
%
%function fp=fpfun(t,x,p)%nxp
%    fp(1,1)=-x(1).*x(3);
%    fp(1,2)=x(2)+x(3)+x(4);
%    fp(1,3)=x(4);
%    fp(2,1)=x(1).*x(3);
%    fp(2,2)=-x(2);
%    fp(2,5)=-x(2);
%    fp(3,2)=-x(3);
%    fp(3,4)=-x(3);
%    fp(3,5)=x(2);
%    fp(4,2)=-x(4);
%    fp(4,3)=-x(4);
%    fp(4,4)=x(3);
%    fp(4,5)=0;
%end

% %SEIRCD
 function dx=ffun(t,x,p)%nx1
     dx(1,1)=-p(1).*exp(-p(3)*t).*x(1).*x(3)./(x(1)+x(2)+x(3)+x(4));
     dx(2,1)=p(1).*exp(-p(3)*t).*x(1).*x(3)./(x(1)+x(2)+x(3)+x(4))-p(4).*x(2);
     dx(3,1)=p(4).*x(2)-p(5).*x(3);
     dx(4,1)=(1-p(2))*p(5).*x(3);
     dx(5,1)=p(4).*x(2);
     dx(6,1)=p(2)*p(5).*x(3);
 end
 
 function fx=fxfun(t,x,p)%nXn
     fx(1,1)=-p(1).*exp(-p(3)*t).*x(3).*(x(2)+x(3)+x(4))./(x(1)+x(2)+x(3)+x(4)).^2;
     fx(1,2)=p(1).*exp(-p(3)*t).*x(3).*x(1)./(x(1)+x(2)+x(3)+x(4)).^2;
     fx(1,3)=-p(1).*exp(-p(3)*t).*x(1).*(x(1)+x(2)+x(4))./(x(1)+x(2)+x(3)+x(4)).^2;
     fx(1,4)=p(1).*exp(-p(3)*t).*x(3).*x(1)./(x(1)+x(2)+x(3)+x(4)).^2;
     fx(2,1)=p(1).*exp(-p(3)*t).*x(3).*(x(2)+x(3)+x(4))./(x(1)+x(2)+x(3)+x(4)).^2;
     fx(2,2)=-p(1).*exp(-p(3)*t).*x(3).*x(1)./(x(1)+x(2)+x(3)+x(4)).^2-p(4);
     fx(2,3)=p(1).*exp(-p(3)*t).*x(1).*(x(1)+x(2)+x(4))./(x(1)+x(2)+x(3)+x(4)).^2;
     fx(2,4)=-p(1).*exp(-p(3)*t).*x(3).*x(1)./(x(1)+x(2)+x(3)+x(4)).^2;
     fx(3,2)=p(4);
     fx(3,3)=-p(5);
     fx(4,3)=(1-p(2))*p(5);
     fx(5,2)=x(4);
     fx(6,3)=p(2)*p(5);
     fx(6,6)=0;
 end
 
 function fp=fpfun(t,x,p)%nxp
     fp(1,1)=-exp(-p(3)*t).*x(1).*x(3)./(x(1)+x(2)+x(3)+x(4));
     fp(1,3)=t*p(1).*exp(-p(3)*t).*x(1).*x(3)./(x(1)+x(2)+x(3)+x(4));
     fp(2,1)=exp(-p(3)*t).*x(1).*x(3)./(x(1)+x(2)+x(3)+x(4));
     fp(2,3)=-t*p(1).*exp(-p(3)*t).*x(1).*x(3)./(x(1)+x(2)+x(3)+x(4));
     fp(2,4)=-x(2);
     fp(3,4)=x(2);
     fp(3,5)=-x(3);
     fp(4,2)=-p(5).*x(3);
     fp(4,5)=(1-p(2)).*x(3);
     fp(5,4)=x(2);
     fp(6,2)=p(5).*x(3);
     fp(6,5)=p(2).*x(3);
 end