%function S= sensitivityanalysis(k,time,xi,S)
clc
clear all

A0 = 0.02090; B0 = A0/3; C0 = 0; D0 = 0; E0 = 0;
data.y0 = [A0;B0;C0;D0;E0];

y0= data.y0;

k(1)=14.402;
k(2)=1.5663;
k(3)=0.29042;

tspan=(0:501/3000:501);
options = odeset('RelTol', 1e-12, 'AbsTol', 1e-12);
[t,z1]=ode45(@(t,y) himmelode(t,y,k),tspan,data.y0,options);

xi=1e-5;
S=[];
I2=[];

for i=1:3
   pad_k=k;
   pad_k(i)=k(i)*(1+xi);
   [t,z2]=ode45(@(t,y) himmelode(t,y,pad_k),tspan,data.y0,options);
    
   mz=mean(z1,1);
   B = repmat(mz,3001,1);
   sens=(z2-z1)/(xi*k(i))*k(i)./B;
   
   S=[S;sens];
   I1=trapz(sens(:,:));
   I2=[I2;trapz(tspan,abs (sens(:,:)))]
   I3=I2'
   
end