lin_evol_opt_fd_fd_Hessian.m

function U_der2 = lin_evol_opt_fd_fd_Hessian(model, model_data, mu)
%function U_der = lin_evol_opt_fd_fd_second_derivative(model, model_data, mu)
%
% function calculating the second derivative solution via PDE and finite difference
% procedure for the parameter mu.
%
% Markus Dihlmann 03.02.2011

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

old_mu = get_mu(model);

model = set_mu(model, mu);

compute_derivative_old = model.compute_derivative_info;

model.compute_derivative_info =1;

%sim_data = detailed_simulation(model, model_data);



U_der_fix=lin_evol_opt_fd_derivative(model, model_data,mu);
U_der2 = cell(length(mu),length(mu));

h=0.0001; %stepsize
for j=1:length(mu)
    for i=1:length(mu)
        %is mu+h inside mu_range?
        if (mu(i)+h>=model.mu_ranges{i}(1))&&(mu(i)+h<=model.mu_ranges{i}(2))
            mu(i) = mu(i) + h;
            model = set_mu(model,mu);
            U_der = lin_evol_opt_fd_derivative(model, model_data,mu);
            mu(i) = mu(i) -h;
            model = set_mu(model,mu);

            U_der2{j,i} = (U_der{j}-U_der_fix{j})./h;
        else
            mu(i) = mu(i) - h;
            model = set_mu(model,mu);
            U_der = lin_evol_opt_fd_derivative(model, model_data,mu);
            mu(i) = mu(i) +h;
            model = set_mu(model,mu);

            U_der2{j,i} = (-U_der{j}+U_der_fix{j})./h;
        end

    end
end

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