function Func = Model(t,y,q1,q2,q3,q4,p1,p2,p3,p4,p5,p6,p7,p8,p9,p10,p11,p12,p13,p14,p15,p16,m1,m2,m3,m4,m5,m6,m7,m8,m9,m10,m11,m12,m13,m14,m15,m16,m17,m18,m19,m20,m21,m22,m23,m24,m25,m26,m27,m28,n0,n1,n2,n3,n4,n5,n6,n7,n8,n9,n10,n11,n12,n13,n14,g1,g2,g3,g4,g5,g6,g7,g8,g9,g10,g11,g12,g13,g14,g15,g16,g17,a,b,c,d,e,f,g,h,i,j,k,l,m,n,o,s,p,dusk,dawn1,dusk1,dusk2,dawn2,day_numb,period);


L=1;
L=0.5*((1+tanh((t-period*floor(t/period))/0.5))-(1+tanh((t-period*floor(t/period)-dusk)/0.5))+(1+tanh((t-period*floor(t/period)-period)/0.5)));   

Func = zeros(21, 1);
%ODEs
% y(1)  LHY mRNA                   
% y(2)  P
% y(3)  GI-ZTL 
% y(4)  
% y(5)  LHY prot
% y(6)  TOC1 mRNA
% y(7)  PRR9 prot
% y(8)  PRR5 (NI) mRNA
% y(9)  PRR5 (NI) prot
% y(10) GI prot
% y(11) TOC1 prot
% y(12) ZTL
% y(13) Y mRNA
% y(14) GI mRNA
% y(15) PRR9 mRNA
% y(16) PRR7 mRNA
% y(17) PRR7 prot
% y(18) TOC1 modif.
% y(19) Y prot.
% y(20) LHY prot modif.
% y(21) 

% 

Func(1) = 1* g1^a/(g1^a+(y(7)+y(17)+y(9))^a)*(n0*L+q1*L*y(2)+n1*y(18)^b/(y(18)^b+g2^b))-y(1)*(m1*L+m2*(1-L));
Func(2) = p7*(1-L)*(1-y(2))-m11*y(2)*L;
Func(3) = p12*L*y(12)*y(10)-p13*y(3)*(1-L)-m21*y(3);
Func(4) = 0;
Func(5) = (p1*L+p2*(1-L))*y(1)-m3*y(5)-p3*y(5)^c/(y(5)^c+g3^c);
Func(6) = 1*(n2*y(19)^d/(y(19)^d+g4^d)+n3)*(g5^e/(g5^e+y(5)^e))-y(6)*m5;
Func(7) = p8*y(15)-(m13*L+m22*(1-L))*y(7);
Func(8) = 1*(n10*y(20)^l/(g12^l+y(20)^l)+n11*y(17)^m/(g13^m+y(17)^m))-m16*y(8);
Func(9) = p10*y(8)-(m17*L+m24*(1-L))*y(9);
Func(10)= p11*y(14)-m19*y(10)-p12*L*y(12)*y(10)+p13*y(3)*(1-L);
Func(11)= p4*y(6)-m8*y(11)-(m6*L+m7*(1-L))*y(11)*(p5*y(12)+y(3));
Func(12)= 1*p14-m20*y(12)-p12*L*y(12)*y(10)+p13*y(3)*(1-L);
Func(13)= 1*(q2*L*y(2)+g7^s/(g7^s+y(11)^s)*g16^g/(g16^g+y(5)^g)*(n5*L+n6*(1-L)))-y(13)*m9;
Func(14)= 1*(q4*L*y(2)+g14^n/(g14^n+y(11)^n)*g15^o/(g15^o+y(5)^o)*n12*L)-y(14)*m18;
Func(15)= 1*(q3*L*y(2)+g8^h/(g8^h+y(11)^h)*(n4+n7*y(5)^i/(y(5)^i+g9^i)))-m12*y(15);
Func(16)= 1*(n8*(y(5)+y(20))^j/(g10^j+(y(5)+y(20))^j)+n9*y(7)^k/(g11^k+y(7)^k))-m14*y(16);
Func(17)= p9*y(16)-y(17)*(m15*L+m23*(1-L));
Func(18)= p15*y(11)^f/(y(11)^f+g6^f)-(m25*L+m26*(1-L))*y(18);
Func(19)= y(13)*p6-m10*y(19);
Func(20)= p3*y(5)^c/(y(5)^c+g3^c)-m4*y(20);
Func(21)= 0;