function [dGdp_fd dGdp_adj] = getFiniteDiffs(m, inp, obj, indsP, indsIC)
% This function is deprecated. However, it is left here because it appears
% that it computes the complex step finite differences, something that does
% not exists otherwise in Kronecker Bio.
opts = [];

% default options
defopts = struct(   'tStop',                {cell(1, length(inp))}, ...
                    'toPlot',                                    1, ...
                    'plotFcn',                            @plotSol, ...
                    'useEvents',             zeros(1, length(inp)), ...
                    'useParams',                             indsP, ...
                    'useICs',                               indsIC, ...
                    'noOpt',                                     0, ...
                    'verbose',                                   0, ...
                    'eval_G_only',                               0, ...
                    'DerivativeCheck',                       'off', ...
                    'tolFiniteDiff',                   [1e-1 1e-3], ...
                    'objTolerance',                           1e-5, ...
                    'tolForward',                      [1e-9 1e-12], ...
                    'tolAdjoint',                      [1e-4 1e-6], ...
                    'tolForwardHess',                  [1e-6 1e-9], ...
                    'tolAdjointHess',                  [1e-6 1e-9], ...
                    'lowerBound',                   [m.p;m.ic]/1e3, ...
                    'upperBound',                   [m.p;m.ic]*1e3, ...
                    'RelLineSrchBnd',                            1, ...
                    'RelLineSrchBndDuration',                  inf, ...
                    'MaxIter',                                 1e3  ...
                    );

opts = mergestruct(defopts, opts);

clear defopts

    % 1a Update the model with the given parameter set
    m           = m.update(m.p, m.ic);

    % 1b Update the objective function
%         obj         = obj(1).update(m);
    % initialize the zero objective function
    obj0 = Gzero(m);
    % augment options structure with default options
    objFields = fieldnames(obj0);
    for i = 1:length(objFields)
        for k = 1:length(obj)
            if ~isfield(obj(k), objFields{i}) || isempty(obj(k).(objFields{i}))
                obj(k).(objFields{i}) = obj0.(objFields{i});
            end
        end
    end

    % 1b update the input
%         inp         = inp(1).update(m);
    % initialize the zero objective function
    inp0 = Uzero(m);
    % augment options structure with default options
    inpFields = fieldnames(inp0);
    for i = 1:length(inpFields)
        for k = 1:length(inp)
            if ~isfield(inp(k), inpFields{i}) || isempty(inp(k).(inpFields{i}))
                inp(k).(inpFields{i}) = inp0.(inpFields{i});
            end
        end
    end
% *******************************************************
%    Forward integration
% *******************************************************
    % simulates and returns objective function value
    disp('integrating forward')
    [G, fwdInfo, deltaTimes] = integrateFwd(m, inp, obj, opts, opts.tolForward);

% *******************************************************
%    Adjoint calculation
% *******************************************************
    % runs if dG/dp is a requested output argument
    disp('integrating adjoint')
    [GradG, adjInfo] = integrateAdj(m, inp, obj, opts, fwdInfo, deltaTimes, opts.tolAdjoint);
    assignin('base', 'backupM', m);
    
% *******************************************************
%    Finite difference integration
% *******************************************************
    p0  = m.p;
    h   = 1e-9;
    for k = 1:length(indsP)
        m.p           = p0;
        m.p(indsP(k))  = m.p(indsP(k)) + sqrt(-1)*h;
        m             = m.update(m.p, m.ic);
        % simulates and returns objective function value
        disp(sprintf('integrating forward for k = %d', k))
        [Gp(k), fwdInfo, deltaTimes] = integrateFwd(m, inp, obj, opts, opts.tolForward);
    end
    m.p = p0;
    ic0 = m.ic;
    for k = 1:length(indsIC)
        m.ic = ic0;
        m.ic(indsIC(k)) = m.ic(indsIC(k)) + sqrt(-1)*h;
        m               = m.update(m.p, m.ic);
        
        % simulates and returns objective function value
        disp(sprintf('integrating forward for k = %d', k))
        [Gp(length(indsP)+k), fwdInfo, deltaTimes] = integrateFwd(m, inp, obj, opts, opts.tolForward);
    end
dGdp_fd = imag(Gp)'/h;
dGdp_adj = GradG';