lin_evol_opt_fd_derivative.m

function U_der = lin_evol_opt_fd_derivative(model, model_data, mu)
%function U_der = lin_evol_opt_fd_derivative(model, model_data, mu)
%
% function calculating the derivative solution via a 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 =0;

sim_data = detailed_simulation(model, model_data);

U_fix=sim_data.U;

h=0.0000001; %stepsize

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_der{i} = (sim_data.U-U_fix)./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_der{i} = (-sim_data.U+U_fix)./h;
    end
    
end

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