lin_evol_opt_fd_second_derivative.m

function U_der2 = lin_evol_opt_fd_second_derivative(model, model_data, mu,h)
%function U_der = lin_evol_opt_fd_second_derivative(model, model_data, mu,h)
%
% 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);



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;
U_der2 = cell(1,length(mu));
if (nargin<4)
    h=0.01; %stepsize
end

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);
        sim_data = detailed_simulation(model, model_data);
        mu(i) = mu(i) -h;
        model = set_mu(model,mu);
        
        U_der2{i} = (sim_data.U_der{i}-U_der_fix{i})./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);
        
        U_der2{i} = (-sim_data.U_der{i}+U_der_fix{i})./h;
    end
    
end

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