lin_evol_opt_fd_Hessian_der.m

function Hessian_der = lin_evol_opt_fd_Hessian_der(model, model_data, mu,h)
%function U_der = lin_evol_opt_fd_Hessian(model, model_data, mu,h)
%
% function caluculating the third derivative of a function, by a
% sensitivity PDE for the first derivative and then finite differences for
% the second derivative of the first derivative :-)
%
% Markus Dihlmann 20.06.2011

if isempty(mu)
    mu = get_mu(model);
end

old_mu = get_mu(model);

if (nargin<4)
    h=0.001; %stepsize
end


compute_derivative_old = model.compute_derivative_info;

model.compute_derivative_info =1;


H= lin_evol_opt_fd_Hessian(model, model_data, mu,h);

n_params=size(H,1);

Hessian_der = cell(size(H,1)^2,n_params);


for k=1:n_params
    mu(k) = mu(k) + h;
    model = set_mu(model, mu);
    Hh = lin_evol_opt_fd_Hessian(model, model_data,mu,h);
    for i=1:n_params
        for j=1:n_params
            Hessian_der{(i-1)*n_params+j,k} = (-H{i,j}+Hh{i,j})./h;
        end
    end
    mu(k) = mu(k) - h;
    model = set_mu(model,mu);
end

model = set_mu(model, old_mu);
model.compute_derivative_info = compute_derivative_old;