//simulations
x0=[-1,1]';dt=linspace(0,20,1001);N=length(dt);t0=dt(1);
a=0.2;b=0.2;c=0.6;p=[a;b;c];disp("p");disp(p);
pourcent=0;//noise in pourcentage

function xdot=fsimu(t,x,p)//The model
     xdot=[p(3)*(x(1)-x(1)^3/3+x(2));-1/p(3)*(x(1)-p(1)+p(2)*x(2))];
endfunction

x=ode(x0,t0,dt,list(fsimu,p));//Simulations
br=rand(ones(N,1));
yb=(1+br.*(0.01*pourcent)).*x(1,:)';//add noise
plot2d(dt,yb,1)//Draw the output

///Program to estimate double integrals
function [Idouble]=Integraledouble(tau,i,dt,functint)//trapezoide method, double integral
    Intsimple=zeros(1,tau);Idouble=0;
    compt=1;
    for j=i-tau:i       
         Intsimple(compt)=inttrap(dt(j-tau:j),functint(j-tau:j));
        compt=compt+1;
    end
    Idouble=inttrap(dt(i-tau:i),Intsimple);
endfunction

//main program. To launch it, do: [A,b,pestim]=EstimationParameters(dt,yb,7); for example 
function [A,b,pestim]=EstimationParameters(dt,y,tau)
   
      //tau is the number of points of integration
      h=dt(2)-dt(1)
    for i=2*tau+2:N 
        //Construction of the vector b 
        b(i)=y(i)-2*y(i-tau)+y(i-2*tau)
        
        //construction of the matrix A   
        A(i,1)=-(tau*h)^2;
        
        vectA2=y.^3;
        A(i,2)=Integraledouble(tau,i,dt,vectA2);
                
        vectA3=y(i-tau:i).^3-y(i-2*tau:i-tau).^3;     
        A3=inttrap(dt(i-tau:i),vectA3)/3;
        A(i,3)=A3;
        
        A(i,4)=Integraledouble(tau,i,dt,y)
         
        vectA5=y(i-tau:i)-y(i-2*tau:i-tau);
        A5=inttrap(dt(i-tau:i),vectA5);
        A(i,5)=A5;
                 
    end
    n1=2*tau+2 ;n2=N-tau
    A1=A(n1:n2,:);[nA,mA]=size(A1);
    disp("conditioning of A1=");disp(cond(A1))
    b1=-b(n1:n2);
    [Q,R]=qr(A1);
    V=Q'*b1;
    R1=zeros(mA);
    for k=1:mA
        R1(k,1:mA)=R(k,1:mA);
    end
    V1=zeros(mA,1);
    V1(1:mA,1)=V(1:mA,1);
    par=zeros(mA,1);
    par=R1\V1;
       
    aestim=par(1); bestim=3*par(2); cestim=par(3);
    
    pestim=[aestim,bestim,cestim]';
    disp('a'); disp(aestim);
    disp('b');disp(bestim);
    disp('c');disp(cestim);
    //Error on the parameter vector
    error1=norm(pestim-p)/norm(p);
    disp("parameter error=");disp(error1);
    //Error on the output
    x=ode(x0,t0,dt,list(fsimu,pestim));
    plot2d(dt,x(1,:)',2)
    error2=norm(y-x(1,:)')/norm(y)   
    disp("output error=");disp(error2);
    
endfunction
