k = ones(1,6)*0.1;
x0 = ones(1,4);

[t,y] = ode45(@dYdT2015, [0 1e3], x0, [], k);
ysol = y(length(y),:);
oldnu = ComputeNu2015(k, ysol);
yss = ysol;

elasticityres = zeros(4,6);
deltay = 0.00002;

for perturb=1:4
    oldy = yss(perturb);
    yss(perturb) = yss(perturb)+deltay;
    [t,y] = ode45(@dYdT2015, [0 1e-6], yss, [], k);
    newysol = y(length(y),:);
    newnu = ComputeNu2015(k, newysol);
    elasticityset = ((newnu-oldnu)/deltay * oldy ./ oldnu);
    for i=1:6
        elasticityres(perturb,i) = elasticityset(i);
    end
    yss(perturb) = oldy;
end