Why is optimization involving ODEs so much slower than with a purely algebraic objective function?

Consider this optimization script, which is a very much boiled down version of something I am currently working on:

% Parameters
substr_0        = 6.21;
substr_in       = 3.171;
k_1             = 0.5;
f_1             = 43;
phi_1           = 0.5;
initialVolume   = 0;

daysTotal       = 7;
hoursTotal      = daysTotal*24;
%% Consumer schedule per hour
onOff = [0 0 0 0 0 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 0 0 0 0 0 0 0 1 1 1 ...
    1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 0 0 0 0 0 0 1 1 1 1 1 1 1 0 0 1 1 1 1 1 ...
    1 1 1 0 0 0 0 0 0 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 0 0 0 0 0 0 0 1 ...
    1 1 1 1 1 1 0 0 0 0 1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 ...
    1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 1 1 1 0]';
% Initial value: Feeding only during "Consumer ON" time
substr_in_time = substr_in/12.*onOff;

%% Setup optimization problem
parameters = [substr_0, k_1, f_1, phi_1, initialVolume];
% Objective function: Standard deviation of storage volume
objFun = @(feedingVector) substr_Model(feedingVector, onOff, parameters);
% Bounds (lb <= x <= ub), i.e., decision variables can only range between 0 and 6.3
lowerBound = zeros(1, hoursTotal);
upperBound = 2*substr_in.*ones(1, hoursTotal);
% Possible solvers: FilterSD, Ipopt, fmincon (local)
optimOpts = optimoptions('fmincon', 'Display', 'iter', 'FunctionTolerance', 1e-2,...
    'MaxIterations', 1e5, 'MaxFunctionEvaluations', 1e5*daysTotal);

%% Solve optimization problem
tic; [feeding_Opt, fval] = fmincon(objFun, substr_in_time, [], [], [], [],...
    lowerBound, upperBound, [], optimOpts); toc

%% Calculate model with optimal values
substr_Model(feeding_Opt, onOff, parameters)

%% Plot results
figure();
plot(1:hoursTotal, interm1, '-o', 'DisplayName', 'Interm production');
hold on; grid on
stairs(0:hoursTotal-1, 5*onOff, '-', 'DisplayName', 'Consumption schedule');
plot(1:hoursTotal, storageVolume/10, '-^', 'MarkerSize', 4,...
    'DisplayName', 'Storage volume [1e-1*m^3]');
stairs(0:hoursTotal-1, feeding_Opt, 'DisplayName', 'Feeding schedule');
xtickVector = 0:12:hoursTotal;
xticks(xtickVector);
xticklabels(xtickVector/24);
xlabel('Time (days)');
legend('Location', 'Best');
set(gca, 'FontSize', 18);
hold off

function storageStd = substr_Model(substr_feeding, onOff, parameters)
    hoursTotal      = length(substr_feeding);
    substr_0        = parameters(1);
    k_1             = parameters(2);
    f_1             = parameters(3);
    phi_1           = parameters(4);
    initialVolume   = parameters(5);

    % Simple for-loop with feeding as vector per hour
    substr = zeros(hoursTotal, 1);
    substr(1) = (substr_0 + substr_feeding(1))*exp(-k_1/24);
    for hour = 2:hoursTotal
       substr(hour) = (substr(hour-1) + substr_feeding(hour))*exp(-k_1/24);
    end
    % Calculate resulting interm1 production and storage volume
    interm1 = phi_1*f_1*substr/24;
    interm1ConsumptionRate = 11.24509;
    storageVolume = initialVolume + cumsum(interm1 - onOff*interm1ConsumptionRate);
    assignin('base', 'substr', substr);
    assignin('base', 'interm1', interm1);
    assignin('base', 'storageVolume', storageVolume);
    storageStd = std(storageVolume);
end

However, If I use an ODE as part of the objective function, optimization takes much longer!

function storageStd = substr_Model(substr_feeding, onOff, parameters)
    hoursTotal      = length(substr_feeding);
    substr_0        = parameters(1);
    k_1             = parameters(2);
    f_1             = parameters(3);
    phi_1           = parameters(4);
    initialVolume   = parameters(5);

    % Simple for-loop with feeding as vector per hour
    substr = zeros(hoursTotal, 1);
    initialValues4Model_0 = [substr_0 + substr_feeding(1)];
    % Solve ODE
    [~, x_sol] = ode15s(@(t, x) substr_Model_diff(t, x, parameters), 0:0.5:1,...
        initialValues4Model_0);
    substr(1) = x_sol(end);
    for hour = 2:hoursTotal
        initialValues4Model = [substr(hour-1) + substr_feeding(hour)];
        [~, x_sol] = ode15s(@(t, x) substr_Model_diff(t, x, parameters), 0:0.5:1,...
            initialValues4Model);
        substr(hour) = x_sol(end);
    end

    % Calculate resulting interm1 production and storage volume
    interm1 = phi_1*f_1*substr/24;
    interm1ConsumptionRate = 11.24509;
    storageVolume = initialVolume + cumsum(interm1 - onOff*interm1ConsumptionRate);
    assignin('base', 'substr', substr);
    assignin('base', 'interm1', interm1);
    assignin('base', 'storageVolume', storageVolume);
    storageStd = std(storageVolume);
end

With the ODE in a separate file:

function dx = substr_Model_diff(t, x, parameters) 
    k_1 = parameters(2);
    % RHS of ODE
    dx	  = -k_1/24*x;
end

(Differential objective code could probably be improved a bit, but this is just for demonstration purposes)

So my question is: Why is optimization involving an ODE so much slower than with a purely algebraic objective function? I understand fmincon uses the obejctive's Jacobian. Is it not able to determine this automatically in the second case? Would passing an analytic Jacobian (created from Symbolic Math Toolbox) help? But then, how to derive it, since the objective contains a for loop?

[dpb]

How could it avoid being "so much slower" when the ODE has to be solved to get the functional value to optimize over?