%s funkcia pre kolonu
function[sys,x0]=kolsf(t,x,u,flag)
f=22.5;
xf=0.23;
g=30; %u=L, G je konstanta a treba pocitat D, pov. L=22.26
l=22.26; %u=G, L je konstanta a treba pocitat D, pov. G=30
%d=7.74; % povodne D
if flag==0,
   sys=[27,0,1,1,0,0];
   xs=[0.0008;0.0030;0.0050;0.0076;0.0111
      0.0156;0.0215;0.0291;0.0387;0.0508
      0.0657;0.0837;0.1050;0.1293;0.1564
      0.1402;0.1419;0.1449;0.1501;0.1592
      0.1749;0.2012;0.2441;0.3106;0.4057
      0.5283;0.6670];   
x0=xs';

elseif flag==1,
h0=1;h=.1;hko=.2;
e=.6;
% pocet etazi a nastrekova etaz
n=25;nf=14;
% koeficienty teoretickej rovnovaznej krivky
ar=.0021567456;cr=2.3066176;er=-0.88823513;
br=0.82112327;dr=-0.39916116;
% vypocet parametrov z celkovej bilancie
d=g-u(1);
%d=u(1)-l;
w=f-d;
%varak - 100% ucinnost
yt(1)=(ar+cr*x(1)+er*x(1)^2)/(1+br*x(1)+dr*x(1)^2);
y(1)=yt(1);
dx(1)=(u(1)+f)/h0*x(2)-g/h0*y(1)-w/h0*x(1);
%1.etaz
%yt(2)=(ar+cr*x(2)+er*x(2)^2)/(1+br*x(2)+dr*x(2)^2);
%y(2)=(1-e)*y(1)+e*yt(2);
%dx(2)=(u(1)+f)/h*(x(3)-x(2))+g/h*(y(1)-y(2))-w/h*x(2);
%1.-nastrekova etaz-1
for i=2:nf
	yt(i)=(ar+cr*x(i)+er*x(i)^2)/(1+br*x(i)+dr*x(i)^2);
	y(i)=(1-e)*y(i-1)+e*yt(i);
	dx(i)=(u(1)+f)/h*(x(i+1)-x(i))+g/h*(y(i-1)-y(i));
end 
% nastrekova etaz
yt(nf+1)=(ar+cr*x(nf+1)+er*x(nf+1)^2)/(1+br*x(nf+1)+dr*x(nf+1)^2);
y(nf+1)=(1-e)*y(nf)+e*yt(nf+1);
dx(nf+1)=f/h*xf+u(1)/h*x(nf+2)-(u(1)+f)/h*x(nf+1)+g/h*(y(nf)-y(nf+1));
%obohacovacia cast
for i=nf+2:n+1
	yt(i)=(ar+cr*x(i)+er*x(i)^2)/(1+br*x(i)+dr*x(i)^2);
	y(i)=(1-e)*y(i-1)+e*yt(i);
	dx(i)=u(1)/h*(x(i+1)-x(i))+g/h*(y(i-1)-y(i));
end 
% posledna etaz
%yt(n+1)=(ar+cr*x(n+1)+er*x(n+1)^2)/(1+br*x(n+1)+dr*x(n+1)^2);
%y(n+1)=(1-e)*y(n)+e*yt(n+1);
%dx(n+1)=u(1)/h*(y(n+1)-x(n+1))+g/h*(y(n)-y(n+1));
% kondenzator
dx(n+2)=g/hko*y(n+1)-(u(1)+d)/hko*x(n+2);

sys=dx;
elseif flag==3,
sys=x(27);
else
sys=[];
end