lin_evol_opt_fd_Hessian.m

function Hessian = lin_evol_opt_fd_Hessian(model, model_data, mu,h)
%function U_der = lin_evol_opt_fd_Hessian(model, model_data, mu,h)
%
% function calculating the second derivative solution via PDE and finite difference
% procedure for the parameter mu.
%
% Markus Dihlmann 20.06.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=sim_data.U_der;
Hessian = cell(length(sim_data.U_der),length(sim_data.U_der));
if (nargin<4)
    h=0.0001; %stepsize
end

for j = 1:length(sim_data.U_der)
    for i=1:length(sim_data.U_der)
        %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);
            sim_data = detailed_simulation(model, model_data);
            mu(i) = mu(i) -h;
            model = set_mu(model,mu);

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

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

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