function dz=getf(t,z)
% Differential equation dz=f(t,z)
% dz(1)=f1(t,z), ... ,dz(n)=fn(t,z) : given equations
% dz(n+1), ... , dz(n+np) : from sensitivity functions

%----- input here -----%
%%%%% write as many sn as the number of parameters
global s1 s2 s3


%----- input here -----%
%%%%% The order of system %%%%%
% n : number of states
% p : number of parameters
n=2;
p=3;


%----- input here -----%
%%%%% Set the variables %%%%%
% x1, x2, ...x(n+np)
% a, b, c, ... p
SYMS  x1 x2 x3 x4 x5 x6 x7 x8 a b c
x(1)=x1;
x(2)=x2;
x(3)=x3;
x(4)=x4;
x(5)=x5;
x(6)=x6;
x(7)=x7;
x(8)=x8;
s(1)=a;
s(2)=b;
s(3)=c;


%----- input here -----%
%%%%% State Equation here %%%%%
% dx(n)=f(x,s)
% state : x1, x2, ... , xn
% parameter : a, b, c, .. p
dx(1)=x2;
dx(2)=-c*sin(x1)-(a+b*cos(x1))*x2;


%%%%% Get df/dx
for i=1:n
    j=1;
    dfdx_=[diff(dx(j),x(i))];
    for j=2:n
        dfdx_=[dfdx_;diff(dx(j),x(i))];
    end    
    dfdx(:,i)=[dfdx_]; 
end

%%%%% Get df/ds
for i=1:p
    j=1;
    dfds_=[diff(dx(j),s(i))];
    for j=2:n
        dfds_=[dfds_;diff(dx(j),s(i))];
    end    
    dfds(:,i)=[dfds_]; 
end

%%%%% construct ODE equations %%%%%
for i=1:p
    j=1;
    S_=[x(n*i+1)];
    for j=2:n
        S_=[S_;x(n+i*j)];
    end    
    S(:,i)=[S_]; 
end

% dS=(df/df)S+(df/ds)
for i=1:p
    for j=1:n
        dx(n*i+j)=subs(dfdx(j,:))*S(:,i)+subs(dfds(j,i));
    end
end


%----- input here -----%
%%%%% get nominal parameter values %%%%%
% a=s1, b=s2, ...
a=s1;
b=s2;
c=s3;
% x1=z(1), x2=z(2),...,x(n+np)=z(n+np)
x1=z(1);
x2=z(2);
x3=z(3);
x4=z(4);
x5=z(5);
x6=z(6);
x7=z(7);
x8=z(8);

dz=subs(dx)';