function [T pakl]=periodas(sistemos_dinamika,apyt_periodas,prad_sal,sprend_tiksl)
% v.003_en_final
% Function periodas() returns vector where the first element is dynamical
% systems period and the second element is period error.
% Functions first argument is dynamical systems equations, second argument
% is approximate period, third argument is initial conditions on limit
% cycle and fourth argument is solution accuracy given the same way like 
% in function ode45  

ck=prad_sal;
T_ap=apyt_periodas;
N= size(ck,2);
% We integrate given system over time span T_ap and save new found point of 
% limit cycle   
[T3 Y3]= ode45(sistemos_dinamika,[0 T_ap],prad_sal,sprend_tiksl);

prad_sal=[Y3(end,:) T_ap];
tikr=false;
interval=T_ap;
clear T3 Y3;
% Now we integrate new system musu_din_sistema2 starting from initial 
% conditions [ck1 T_ap] as long as a solution saturate. 

while tikr~=true

[T1 Y1]= ode45(@musu_din_sistema2,[0 interval],prad_sal,sprend_tiksl); 
    if abs(Y1(end,N+1)-Y1(end-1,N+1))<min([odeget(sprend_tiksl,'Reltol') odeget(sprend_tiksl,'AbsTol')])   
    per=Y1(end,N+1);
    tikr=true;    
    else      
    prad_sal=Y1(end,:);
    interval=T_ap;
    end;  
    
end;

function dydt=musu_din_sistema2(t,y)   
   dydt(1:N,1)=-2*((y(1:N)'-ck)*sistemos_dinamika(t,y(1:N)))*sistemos_dinamika(t,y(1:N));
   dydt(N+1,1)=-2*((y(1:N)'-ck)*sistemos_dinamika(t,y(1:N)));
end

T=per;
pakl=abs(Y1(end,N+1)-Y1(end-1,N+1));


end