delta=10^(-4);
iter=5;        %number of iterations 
k=zeros(iter+1,np);
fac=0.5;
k1=g1*fac;
k2=g2*fac;
k3=g3*fac;
%[ms(g,:);mt(g,:);mu(g,:);mv(g,:)];

%%%%%%%%%%%%%%%%%%%%
k(1,1)=k1;
k(1,2)=k2;
k(1,3)=k3;
%%%%%%%%%%%%%%%%%%%%
uobs=1:obs+1;
%[1 2 4 5 7 9 10 11 12 13 14 15 16];
nobs=length(uobs)-1;
%%%%%%%%%%%%%%%%%%%%
ms=zeros(iter+1,nobs+1);
mt=zeros(iter+1,nobs+1);
mu=zeros(iter+1,nobs+1);
mv=zeros(iter+1,nobs+1);
mes=zeros(nobs+1,1);
met=zeros(nobs+1,1);
meu=zeros(nobs+1,1);
mev=zeros(nobs+1,1);
b=zeros(ne*nobs,1);
bb=zeros(ne*nobs,1);
d=zeros(ne*(nobs+1),1);
dd=zeros(ne*(nobs+1),1);
ss=zeros(iter,nobs+1+np,nobs);
tt=zeros(iter,nobs+1+np,nobs);
uu=zeros(iter,nobs+1+np,nobs);
vv=zeros(iter,nobs+1+np,nobs);
A=zeros(ne*nobs,ne*(nobs+1)+np);
timeseries1=zeros(iter,nobs*res);
timeseries2=zeros(iter,nobs*res);
timeseries3=zeros(iter,nobs*res);
timeseries4=zeros(iter,nobs*res);
tp=zeros(1,nobs+1);
ix=zeros(1,nobs+1);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
s=stuv(:,1);
t=stuv(:,2);
u=stuv(:,3);
v=stuv(:,4);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

for i=1:nobs+1
    tp(i)=t0+h*res*(uobs(i)-1);
    ix(i)=1+res*(uobs(i)-1);
    mes(i)=s(uobs(i));
    met(i)=t(uobs(i));
    meu(i)=u(uobs(i));
    mev(i)=v(uobs(i));
    ms(1,i)=mes(i);
    mt(1,i)=met(i);
    mu(1,i)=meu(i);
    mv(1,i)=mev(i);
end;


ind=1;
trh=10^(-7);
g=0;
while ind>0,
    g=g+1
    for j=1:nobs
        [tim,y] = ode45(@zeq,[tp(j):h:tp(j+1)],[ms(g,j) mt(g,j) mu(g,j) mv(g,j)],[],k(g,1),k(g,2),k(g,3));
        timeseries1(g,ix(j):ix(j+1)-1)=y(1:end-1,1)';
        timeseries2(g,ix(j):ix(j+1)-1)=y(1:end-1,2)';
        timeseries3(g,ix(j):ix(j+1)-1)=y(1:end-1,3)';
        timeseries4(g,ix(j):ix(j+1)-1)=y(1:end-1,4)';
        ss(g,1,j)=y(end,1);
        tt(g,1,j)=y(end,2);
        uu(g,1,j)=y(end,3);
        vv(g,1,j)=y(end,4);
        b(j)=ms(g,j+1)-ss(g,1,j);
        b(j+nobs)=mt(g,j+1)-tt(g,1,j);
        b(j+2*nobs)=mu(g,j+1)-uu(g,1,j);
        b(j+3*nobs)=mv(g,j+1)-vv(g,1,j);
        d(j)=mes(j+1)-ss(g,1,j);
        d(j+nobs)=met(j+1)-tt(g,1,j);
        d(j+2*nobs)=meu(j+1)-uu(g,1,j);
        d(j+3*nobs)=mev(j+1)-vv(g,1,j);
    end;
    %k1%%%%%%%%%%%%%%
    for j=1:nobs
        [tim,y] = ode45(@zeq,[tp(j):h:tp(j+1)],[ms(g,j) mt(g,j) mu(g,j) mv(g,j)],[],k(g,1)+delta,k(g,2),k(g,3));
        ss(g,2,j)=y(end,1);
        tt(g,2,j)=y(end,2);
        uu(g,2,j)=y(end,3);
        vv(g,2,j)=y(end,4);
        A(j,1)=(ss(g,2,j)-ss(g,1,j))/delta;
        A(nobs+j,1)=(tt(g,2,j)-tt(g,1,j))/delta;
        A(2*nobs+j,1)=(uu(g,2,j)-uu(g,1,j))/delta;
        A(3*nobs+j,1)=(vv(g,2,j)-vv(g,1,j))/delta;
    end;
    %k2%%%%%%%%%%%%%%
    for j=1:nobs
        [tim,y] = ode45(@zeq,[tp(j):h:tp(j+1)],[ms(g,j) mt(g,j) mu(g,j) mv(g,j)],[],k(g,1),k(g,2)+delta,k(g,3));
        ss(g,3,j)=y(end,1);
        tt(g,3,j)=y(end,2);
        uu(g,3,j)=y(end,3);
        vv(g,3,j)=y(end,4);
        A(j,2)=(ss(g,3,j)-ss(g,1,j))/delta;
        A(nobs+j,2)=(tt(g,3,j)-tt(g,1,j))/delta;
        A(2*nobs+j,2)=(uu(g,3,j)-uu(g,1,j))/delta;
        A(3*nobs+j,2)=(vv(g,3,j)-vv(g,1,j))/delta;
    end;
    %k3%%%%%%%%%%%%%
    for j=1:nobs
        [tim,y] = ode45(@zeq,[tp(j):h:tp(j+1)],[ms(g,j) mt(g,j) mu(g,j) mv(g,j)],[],k(g,1),k(g,2),k(g,3)+delta);
        ss(g,4,j)=y(end,1);
        tt(g,4,j)=y(end,2);
        uu(g,4,j)=y(end,3);
        vv(g,4,j)=y(end,4);
        A(j,3)=(ss(g,4,j)-ss(g,1,j))/delta;
        A(nobs+j,3)=(tt(g,4,j)-tt(g,1,j))/delta;
        A(2*nobs+j,3)=(uu(g,4,j)-uu(g,1,j))/delta;
        A(3*nobs+j,3)=(vv(g,4,j)-vv(g,1,j))/delta;
    end;
    %ss%%%%%%%%%%%%%
    for j=1:nobs
        [tim,y] = ode45(@zeq,[tp(j):h:tp(j+1)],[ms(g,j)+delta mt(g,j) mu(g,j) mv(g,j)],[],k(g,1),k(g,2),k(g,3));
        ss(g,1+np+j,j)=y(end,1);
        tt(g,1+np+j,j)=y(end,2);
        uu(g,1+np+j,j)=y(end,3);
        vv(g,1+np+j,j)=y(end,4);
        A(j,np+j)=(ss(g,1+np+j,j)-ss(g,1,j))/delta;
        A(nobs+j,np+j)=(tt(g,1+np+j,j)-tt(g,1,j))/delta;
        A(2*nobs+j,np+j)=(uu(g,1+np+j,j)-uu(g,1,j))/delta;
        A(3*nobs+j,np+j)=(vv(g,1+np+j,j)-vv(g,1,j))/delta;
        A(j,1+np+j)=-1;
    end;
    %tt%%%%%%%%%%%%%
    for j=1:nobs
        [tim,y] = ode45(@zeq,[tp(j):h:tp(j+1)],[ms(g,j) mt(g,j)+delta mu(g,j) mv(g,j)],[],k(g,1),k(g,2),k(g,3));
        ss(g,1+np+j,j)=y(end,1);
        tt(g,1+np+j,j)=y(end,2);
        uu(g,1+np+j,j)=y(end,3);
        vv(g,1+np+j,j)=y(end,4);
        A(j,nobs+1+np+j)=(ss(g,1+np+j,j)-ss(g,1,j))/delta;
        A(nobs+j,nobs+1+np+j)=(tt(g,1+np+j,j)-tt(g,1,j))/delta;
        A(2*nobs+j,nobs+1+np+j)=(uu(g,1+np+j,j)-uu(g,1,j))/delta;
        A(3*nobs+j,nobs+1+np+j)=(vv(g,1+np+j,j)-vv(g,1,j))/delta;
        A(nobs+j,nobs+2+np+j)=-1;
    end;   
    %uu%%%%%%%%%%%%%
    for j=1:nobs
        [tim,y] = ode45(@zeq,[tp(j):h:tp(j+1)],[ms(g,j) mt(g,j) mu(g,j)+delta mv(g,j)],[],k(g,1),k(g,2),k(g,3));
        ss(g,1+np+j,j)=y(end,1);
        tt(g,1+np+j,j)=y(end,2);
        uu(g,1+np+j,j)=y(end,3);
        vv(g,1+np+j,j)=y(end,4);
        A(j,2*(nobs+1)+np+j)=(ss(g,1+np+j,j)-ss(g,1,j))/delta;
        A(nobs+j,2*(nobs+1)+np+j)=(tt(g,1+np+j,j)-tt(g,1,j))/delta;
        A(2*nobs+j,2*(nobs+1)+np+j)=(uu(g,1+np+j,j)-uu(g,1,j))/delta;
        A(3*nobs+j,2*(nobs+1)+np+j)=(vv(g,1+np+j,j)-vv(g,1,j))/delta;
        A(2*nobs+j,2*(nobs+1)+1+np+j)=-1;
    end;  
    %vv%%%%%%%%%%%%%
    for j=1:nobs
        [tim,y] = ode45(@zeq,[tp(j):h:tp(j+1)],[ms(g,j) mt(g,j) mu(g,j) mv(g,j)+delta],[],k(g,1),k(g,2),k(g,3));
        ss(g,1+np+j,j)=y(end,1);
        tt(g,1+np+j,j)=y(end,2);
        uu(g,1+np+j,j)=y(end,3);
        vv(g,1+np+j,j)=y(end,4);
        A(j,3*(nobs+1)+np+j)=(ss(g,1+np+j,j)-ss(g,1,j))/delta;
        A(nobs+j,3*(nobs+1)+np+j)=(tt(g,1+np+j,j)-tt(g,1,j))/delta;
        A(2*nobs+j,3*(nobs+1)+np+j)=(uu(g,1+np+j,j)-uu(g,1,j))/delta;
        A(3*nobs+j,3*(nobs+1)+np+j)=(vv(g,1+np+j,j)-vv(g,1,j))/delta;
        A(3*nobs+j,3*(nobs+1)+1+np+j)=-1;
    end;  
    %%%%%%%%%%%%%%%%%
    C=A;
    for w=1:ne
        for i=1:nobs
            C((w-1)*nobs+i,(w-1)*(nobs+1)+i+np+1)=0;
        end;
    end;
    
    C=[zeros(1,np),1,zeros(1,4*(nobs+1)-1);C(1:nobs,:);
       zeros(1,np+nobs+1),1,zeros(1,3*(nobs+1)-1);C(nobs+1:2*nobs,:);
       zeros(1,np+2*(nobs+1)),1,zeros(1,2*(nobs+1)-1);C(2*nobs+1:3*nobs,:);
       zeros(1,np+3*(nobs+1)),1,zeros(1,(nobs+1)-1);C(3*nobs+1:4*nobs,:)];
    d=[mes(1)-ms(g,1);d(1:nobs);
       met(1)-mt(g,1);d(nobs+1:2*nobs);
       meu(1)-mu(g,1);d(2*nobs+1:3*nobs);
       mev(1)-mv(g,1);d(3*nobs+1:4*nobs)];
    E=[A;C];
    f=[b;d];
    %%%%%%%%%%%%%%%%%%%%%
    %be=b(1:nobs);
    %Aa=A(1:obs,:);
    %%%%%%%%%%%%%%%%%%%%%
    sol = lsqlin(C,d,[],[],A,b);
    o=1;
    p=0;
    while o>0,
    
    for r=1:np
        k(g+1,r)=k(g,r)+2^(-p)*sol(r);
    end;
    
    ms(g+1,:)=ms(g,:)+2^(-p)*sol(np+1:np+1+nobs)';
    mt(g+1,:)=mt(g,:)+2^(-p)*sol(np+2+nobs:2*(nobs+1)+np)';
    mu(g+1,:)=mu(g,:)+2^(-p)*sol(2*(nobs+1)+np+1:3*(nobs+1)+np)';
    mv(g+1,:)=mv(g,:)+2^(-p)*sol(3*(nobs+1)+np+1:4*(nobs+1)+np)';
    
    for j=1:nobs
        [tim,y] = ode45(@zeq,[tp(j):h:tp(j+1)],[ms(g+1,j) mt(g+1,j) mu(g+1,j) mv(g+1,j)],[],k(g+1,1),k(g+1,2),k(g+1,3));
        ss(g,1,j)=y(end,1);
        tt(g,1,j)=y(end,2);
        uu(g,1,j)=y(end,3);
        vv(g,1,j)=y(end,4);
        bb(j)=ms(g+1,j+1)-ss(g,1,j);
        bb(j+nobs)=mt(g+1,j+1)-tt(g,1,j);
        bb(j+2*nobs)=mu(g+1,j+1)-uu(g,1,j);
        bb(j+3*nobs)=mv(g+1,j+1)-vv(g,1,j);
        dd(j)=mes(j+1)-ss(g,1,j);
        dd(j+nobs)=met(j+1)-tt(g,1,j);
        dd(j+2*nobs)=meu(j+1)-uu(g,1,j);
        dd(j+3*nobs)=mev(j+1)-vv(g,1,j);
    end;
    dd=[mes(1)-ms(g+1,1);dd(1:nobs);
       met(1)-mt(g+1,1);dd(nobs+1:2*nobs);
       meu(1)-mu(g+1,1);dd(2*nobs+1:3*nobs);
       mev(1)-mv(g+1,1);dd(3*nobs+1:4*nobs)];
    ff=[bb;dd];
    
    if or(norm(E\ff,2)<norm(E\f,2),p>=10)
        %(and(norm(bb,2)<norm(b,2),norm(dd,2)<norm(d,2)),p>=10)
       o=0;
       [norm(b,2) norm(d,2)]
    end;
    p=p+1;
    end;
    
    if or(norm(b,2)<trh,g>=iter)
        ind=0;
    end;
end;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


ixx=zeros(nobs+1,1);
for i=1:nobs+1
    ixx(i)=1+res*(uobs(i)-1);
end;

%for i=1:obs+1
%    mod(i,1)=timeseries1(g,min(1+res*(i-1),step));
%    mod(i,2)=timeseries2(g,min(1+res*(i-1),step));
%    mod(i,3)=timeseries3(g,min(1+res*(i-1),step));
%    mod(i,4)=timeseries4(g,min(1+res*(i-1),step));
%end;
%mod;
ido=1:2100;

plot(ido/210,timeseries1(g,:),'k-','linewidth',2);
hold on;
%plot(ido/210,x(1:2100,1),'k-');
hold on;
plot(ixx/210,mes,'bo','linewidth',2);
hold on;
plot(ido/210,timeseries2(g,:),'k:','linewidth',2);
hold on;
%plot(ido/210,x(1:2100,2),'k:');
hold on;
plot(ixx/210,met,'b+','linewidth',2);
hold on;
plot(ido/210,timeseries3(g,:),'k-.','linewidth',2);
hold on;
%plot(ido/210,x(1:2100,3),'k-.');
hold on;
plot(ixx/210,meu,'b*','linewidth',2);
hold on;
plot(ido/210,timeseries4(g,:),'k--','linewidth',2);
hold on;
%plot(ido/210,x(1:2100,4),'b--');
hold on;
plot(ixx/210,mev,'bv','linewidth',2);
title('Enzyme-kinetic reaction');
ylabel('concentration of proteins [uM]'), xlabel('time [min]');
legend('-','Substrate',':','Enzyme','-.','Complex','--','Product');
axis([0 10 0 12]);


%x1%
%plot(ido/210,timeseries3(1,:),'k.','linewidth',2);
%hold on;
%plot(ixx/210,meu,'b*','linewidth',2);


%plot(ido/210,ideig,'r.','linewidth',2);
%hold on;
%plot(ixx/210,meu,'b*','linewidth',2);


%%%%%%%%%%%%%%%%%%%%%%%%%%
k(1:g+1,:)
