function [stock] = stock_ode(flow,timespan_flow,...
    pdf,timespan_pdf,pdf_scale_factor,timesteps_pdf_scale_factor,timespan_out)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%
%%% Reconstructing stock from the ODE
%%%
%%% ds/dt = POM - Integral(POM)
%%%
%%% The Integral is taken back in time, with POM weighted with a given pdf (POM -> WG
%%% through residence time, etc), and with the pdf itself varying with time:
%%% pdf = rho(t-t';t'); all the entries in the pdf matrix are assumed to
%%% have been foramted to the same timespan, while t' covers the entire 
%%% timespan_out
%%%
%%% Dmitry Yumashev, 22/12/2019 
%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%% working out the integral term

[integral, flow_extrap] = probability_integral_advanced(flow,timespan_flow,...
    pdf,timespan_pdf,pdf_scale_factor,timesteps_pdf_scale_factor,timespan_out);


%%% make sure the time dimensions are correct to allow for matrix
%%% operations
if size(timespan_out,1) > 1 && size(timespan_out,2) == 1 % if the vector variable is a column
    timespan_out = transpose(timespan_out); % make it a row
end


%%% general settings

total_timesteps_out = numel(timespan_out);
stock = zeros(1,total_timesteps_out); % make sure it's a row


%%% now solve the ode for the stock

stock(1) = 0; % initial condition: always start with zero stock

for i=2:total_timesteps_out % NOTE the starting index

    stock(i) = stock(i-1) + (flow_extrap(i-1) - integral(i-1)) * (timespan_out(i) - timespan_out(i-1));
    
end
    
end