%g1 = 0.18;
%g2 = 0.02;  %real parameters
%g3 = 0.23;
iter=g;     %use of 'g' will chenge in this program
g1=k(iter,1);
g2=k(iter,2);   %obtained parameters
g3=k(iter,3);
np=3;       %number of parameters
ne=4;       %number of equations
sdev=0.2;
x0=12;
y0=12;
z0=0;
w0=0;
t0=0;
t1=10;
step=2100;
obs=15;     %observations minus one
h=(t1-t0)/step;
res=step/obs;
delta=10^(-4);
% === ODE's method =========
[t,x] = ode45(@zeq,[t0:h:t1+h],[x0 y0 z0 w0],[],g1,g2,g3);
f=[-g1*x(:,1).*x(:,2) + g2*x(:,3), -g1*x(:,1).*x(:,2) + g2*x(:,3) + g3*x(:,3), g1*x(:,1).*x(:,2) - g2*x(:,3) - g3*x(:,3), g3*x(:,3)];
[t,a] = ode45(@zeq,[t0:h:t1+h],[x0+delta y0 z0 w0],[],g1,g2,g3);
ka=[-g1*a(:,1).*a(:,2) + g2*a(:,3), -g1*a(:,1).*a(:,2) + g2*a(:,3) + g3*a(:,3), g1*a(:,1).*a(:,2) - g2*a(:,3) - g3*a(:,3), g3*a(:,3)];
[t,b] = ode45(@zeq,[t0:h:t1+h],[x0 y0+delta z0 w0],[],g1,g2,g3);
kb=[-g1*b(:,1).*b(:,2) + g2*b(:,3), -g1*b(:,1).*b(:,2) + g2*b(:,3) + g3*b(:,3), g1*b(:,1).*b(:,2) - g2*b(:,3) - g3*b(:,3), g3*b(:,3)];
[t,c] = ode45(@zeq,[t0:h:t1+h],[x0 y0 z0+delta w0],[],g1,g2,g3);
kc=[-g1*c(:,1).*c(:,2) + g2*c(:,3), -g1*c(:,1).*c(:,2) + g2*c(:,3) + g3*c(:,3), g1*c(:,1).*c(:,2) - g2*c(:,3) - g3*c(:,3), g3*c(:,3)];
[t,d] = ode45(@zeq,[t0:h:t1+h],[x0 y0 z0 w0+delta],[],g1,g2,g3);
kd=[-g1*d(:,1).*d(:,2) + g2*d(:,3), -g1*d(:,1).*d(:,2) + g2*d(:,3) + g3*d(:,3), g1*d(:,1).*d(:,2) - g2*d(:,3) - g3*d(:,3), g3*d(:,3)];
[t,e] = ode45(@zeq,[t0:h:t1+h],[x0 y0 z0 w0],[],g1+delta,g2,g3);
ke=[-g1*e(:,1).*e(:,2) + g2*e(:,3), -g1*e(:,1).*e(:,2) + g2*e(:,3) + g3*e(:,3), g1*e(:,1).*e(:,2) - g2*e(:,3) - g3*e(:,3), g3*e(:,3)];
[t,f] = ode45(@zeq,[t0:h:t1+h],[x0 y0 z0 w0],[],g1,g2+delta,g3);
kf=[-g1*f(:,1).*f(:,2) + g2*f(:,3), -g1*f(:,1).*f(:,2) + g2*f(:,3) + g3*f(:,3), g1*f(:,1).*f(:,2) - g2*f(:,3) - g3*f(:,3), g3*f(:,3)];
[t,g] = ode45(@zeq,[t0:h:t1+h],[x0 y0 z0 w0],[],g1,g2,g3+delta);
kg=[-g1*g(:,1).*g(:,2) + g2*g(:,3), -g1*g(:,1).*g(:,2) + g2*g(:,3) + g3*g(:,3), g1*g(:,1).*g(:,2) - g2*g(:,3) - g3*g(:,3), g3*g(:,3)];
aa=(a-x)/delta;
bb=(b-x)/delta;
cc=(c-x)/delta;
dd=(d-x)/delta;
ee=(e-x)/delta;
ff=(f-x)/delta;
gg=(g-x)/delta;
kaa=(ka-f)/delta;
kbb=(kb-f)/delta;
kcc=(kc-f)/delta;
kdd=(kd-f)/delta;
kee=(ke-f)/delta;
kff=(kf-f)/delta;
kgg=(kg-f)/delta;
%%%%%%%%%%%%%%%%%%%%
tim=1:res:step+1;
sum=zeros(ne+np);
%%%%%%%%%%
%tim=[1 1 1 1 1 75 76 89 89 89 697 697 796 796 796 1118];
%%%%%%%%%

for j=1:obs+1
    for i=1:ne
        sum=sum+[aa(tim(j),i);bb(tim(j),i);cc(tim(j),i);dd(tim(j),i);ee(tim(j),i);ff(tim(j),i);gg(tim(j),i)]*[aa(tim(j),i);bb(tim(j),i);cc(tim(j),i);dd(tim(j),i);ee(tim(j),i);ff(tim(j),i);gg(tim(j),i)]';
    end;
end;
eF=sum;
eeF=eF;
covar=inv(sum);
det(eF)
%score0=covar(5,5)/g1^2+covar(6,6)/g2^2+covar(7,7)/g3^2
%%%%%%%%%%%%%%%%%%%%
%ests=ms(g,:);
%estt=mt(g,:);
%estu=mu(g,:);
%estv=mv(g,:);
%ep=stuv-[ests' estt' estu' estv'];
%ep=stuv - real;
%b0=zeros(ne+np,1);
%for j=1:obs+1
%    for i=1:ne
%        b0=b0+ep(j,i)*[aa(tim(j),i);bb(tim(j),i);cc(tim(j),i);dd(tim(j),i);ee(tim(j),i);ff(tim(j),i);gg(tim(j),i)];
%    end;
%end;
%dteta=covar*b0;
%dk=dteta(5:7)./[g1;g2;g3];
%norm(dk,2)
%noi=zeros(1:16,1);
%for i=1:16
%    noi(i)=norm(ep(i,:),2);
%end;
%noi

