MATLAB: How to handle discrete, non-periodic right-hand side of ODE in Matlab

differential equationsdiscontinuityMATLABnumerical integrationode15ssimulink

I am trying to migrate an existing Simulink simulation to pure Matlab, which is necessary for several reasons. It is a stiff system of 10 – 50 ordinary, first-order, non-linear differential equations and a couple of algebraic equations, formulated as a level-2 C-Code S-function.
A simplified version of the corresponding Simulink flowsheet (strangely enough, I cannot upload slx files here):
Which produces the following plot:
Consider this MWE for illustration purposes (Matlab version of the Simulink sheet):
% simpleStorageTankModel_Matlab Test of Matlab Solver for Simulink problem
contOutlet = @(t) 0.03*t.^2 - 0.05*t + 0.01;
t1 = 0:1e-2:12;
stairs(t1, periodicDiscreteInlet(t1), 'b'); hold on
% Formulate ODE
solverOptions = odeset('MaxStep', 1e-2);
tic; fillingLevel = ode15s(@(t,y) periodicDiscreteInlet(t) - contOutlet(t),...
    t1, 0.5, solverOptions); toc
plot(t1, 0.03*t1.^2 - 0.05*t1 + 0.01, 'k');
plot(fillingLevel.x, fillingLevel.y, 'r');
legend({'Inlet', 'Outlet', 'Filling level'}, 'Location', 'NorthWest');
function perDisInp = periodicDiscreteInlet(t)
    amplitude = 4;
    phaseDelay = 0;
    pulseWidth = 0.5;
    t_cyc = t - floor(t);
    if numel(t) > 1
        numberOfTimePoints = length(t_cyc);
        perDisInp = zeros(numberOfTimePoints, 1);
        for timeIndex = 1:numberOfTimePoints
            if (t_cyc(timeIndex) >= phaseDelay) && (t_cyc(timeIndex) <...
                1 - pulseWidth + phaseDelay)
                perDisInp(timeIndex) = amplitude;
            end
        end
    else
        if (t_cyc >= phaseDelay) && (t_cyc < 1 - pulseWidth + phaseDelay)
            perDisInp = amplitude;
        else
            perDisInp = 0;
        end
    end
end
An annoying problem is that, if do not set the MaxStep option to something below, say, 1e-2, then the solution beyond t > 6 becomes inaccurate, as if the solver is missing some of the square waves. This seems to have occured to others before .
However, my main problem is with ode15s ignoring the time vector argument I am providing, probably due to the right-hand side of the ODE ( periodicDiscreteInlet) being discrete, i.e., not continuously differentiable (in my simulation, the input is actually a non-periodic square wave). Thus, although t1 is a 1×1201 double vector, I get a 1×1224 double solution. I know I could write a function to extract the values I need, but this would add a significant performance penalty. In my case, Simulink takes just about a tenth of the simulation time of pure Matlab!
Is there a built-in function to "round off" the square-wave edges, making the input continuous (I don't mean as round as via Fourier transform)? Or a stiff Matlab solver that can handle discrete right-hand sides? Why can Simulink handle this problem with such ease, but Matlab can't? Any help is appreciated.

Best Answer

Since this was inspired by the valuable comment of Walter Roberson to my original answer, I am posting this solution as a separate answer.
So it turns out the for-loop solution actually performs 6-8x faster than interp1 inside the ODE RHS function! Here's a plot of both solutions:
For-loop vs. interp1
Code below.
% simpleStorageTankModel_Matlab Test of Matlab Solver for Simulink problem
clc; clearvars
p1 = [0.03 -0.05 0.01]; % contOutlet = @(t) 0.03*t.^2 - 0.05*t + 0.01;
t1 = 0:0.5:12;
pdi = periodicDiscreteInlet(t1);
discreteInput = [t1', pdi, 0.5*ones(25,1)];
%% Solve ODE in a loop without interp1
fillingLevel1 = zeros(length(t1), 1);
initialFillingLevel = 0.5;
fillingLevel1(1) = initialFillingLevel;
tic
for timeInd = 1:length(t1)-1
    odeRhs1 = @(t, y) periodicDiscreteInlet(t1(timeInd)) - polyval(p1, t);
    [~, solVec] = ode15s(odeRhs1, linspace(t1(timeInd), t1(timeInd+1), 3), ...
        initialFillingLevel);
    fillingLevel1(timeInd+1) = solVec(end); % Returns 3 values
    initialFillingLevel = solVec(end);
end
toc
%% Solve ODE with interp1
solverOptions = odeset('MaxStep', 1e-2);
modelOdeFun = @(t, y) odeRhs2(t, y, discreteInput, p1);
tic; [~, fillingLevel2] = ode15s(modelOdeFun, t1, 0.5, solverOptions); toc
%% Plot results: V1 as loop, V2 with interpolation
set(groot, 'DefaultLineLineWidth', 1.25);
set(groot, 'DefaultStairLineWidth', 1.25);
figure();
stairs(discreteInput(:,1), discreteInput(:,2), 'b');
hold on; grid on
plot(t1, polyval(p1, t1), 'k');
plot(t1, fillingLevel1, 'r--');
plot(t1, fillingLevel2);
legend({'Inlet', 'Outlet', 'Filling level-loop', 'Filling level-interp1'}, ...
    'Location', 'NorthWest');
title('RHS with discrete jumps, loop vs. interpolation');
set(gca, 'FontSize', 12);
%% Separate functions
function dy = odeRhs2(t, y, discreteInput, p1)
    % Interpolation as specified in Matlab documentation
    % https://mathworks.com/help/matlab/ref/ode45.html#bu3l43b 
    discreteInput_t = interp1(discreteInput(:,1), discreteInput(:, 2), t, 'previous', 0);
    
    dy = discreteInput_t - polyval(p1, t);
end
function perDisInp = periodicDiscreteInlet(t)
    amplitude = 4;
    phaseDelay = 0;
    pulseWidth = 0.5;
    t_cyc = t - floor(t);
    if numel(t) > 1
        numberOfTimePoints = length(t_cyc);
        perDisInp = zeros(numberOfTimePoints, 1);
        for timeIndex = 1:numberOfTimePoints
            if (t_cyc(timeIndex) >= phaseDelay) && (t_cyc(timeIndex) <...
                    1 - pulseWidth + phaseDelay)
                perDisInp(timeIndex) = amplitude;
            end
        end
    else
        if (t_cyc >= phaseDelay) && (t_cyc < 1 - pulseWidth + phaseDelay)
            perDisInp = amplitude;
        else
            perDisInp = 0;
        end
    end
end
Related Question