Here is an example implementation of the Probability Distributed Model.

function [solution, hydrograph, obj] = simulate(obj, P, t_max, nt)
      
      % Calculate length of each time step
      dt = t_max / nt;
      
      % If precipitation is specified with a single value set precipitation
      % rate to be the same for all time steps
      if length(P) == 1
        P = P * ones(1, nt);
      end
      
      % Initialise structure to save store values (an extra entry is
      % included to store the initial values for each store)
      
      solution.C = zeros(1, nt+1);  % soil moisture storage capacity
      solution.S = zeros(1, nt+1);  % soil moisture storage value
      solution.Sf = zeros(1, nt+1); % fast storage value
      solution.Ss = zeros(1, nt+1); % slow storage value
      
      % Update stores in every time step
      for i = 1:nt+1
        
        % Calculate the soil mosisture value based on the current soil
        % moisture capacity

        if obj.C < obj.par.c_max
          S = obj.par.S_max * (1 - (1 - obj.C / obj.par.c_max) .^ (obj.par.b + 1));
        else
          S = obj.par.S_max;
        end
        
        % Save store values (for i=1 initial value is saved)
        solution.C(i) = obj.C;
        solution.S(i) = S;
        solution.Sf(i) = obj.Sf;
        solution.Ss(i) = obj.Ss;
        
        % If the last step was reached end the simulation
        if i == nt + 1
          break
        end
        
        % Calculate the drainage function
        drainage = max(0, (S - obj.par.st) / obj.par.kg);
        
        % Calculate Pi (source term in ODE for the soil moisture storage)
        Pi = P(i) - drainage;
        
        % Save old value of c (for runoff calculation)
        C_old = obj.C;
      
        % Update value of soil moisture capacity, c
        obj.C = min(obj.C + Pi * dt, obj.par.c_max);
        
        % Calculate the volume of the surface runoff from the water balance
        runoff = Pi * dt - obj.par.S_max * ...
          ((1 - C_old / obj.par.c_max) ^ (obj.par.b + 1) - ...
          (1 - obj.C / obj.par.c_max) ^ (obj.par.b + 1));
        
        % Calculare the rate of runoff generation
        runoff = runoff / dt;
        
        % Calculate slow and fast flow
        flow_s = exp(obj.Ss / obj.par.ks);
        flow_f = obj.Sf / obj.par.kf;
        
        % Update values of the slow and fast storage
        obj.Sf = obj.Sf + (runoff - flow_f) * dt;
        obj.Ss = obj.Ss + (drainage - flow_s) * dt;
      end
      
      % For each time step calculate the corresponding time (t),
      % fast flow (Qf), slow flow (Qs) and total flow (Q_total)
      hydrograph.t = linspace(0, t_max, nt+1);
      hydrograph.Qf = solution.Sf / obj.par.kf;
      hydrograph.Qs = exp(solution.Ss / obj.par.ks);
      hydrograph.Q_total = hydrograph.Qs + hydrograph.Qf;
    end