Level 2 S function

I want to integrate a system of ten ordinary differential equations (ODE's) numerically. Usually, this could be done by for example the ode45-solver. In my case, to compute the derivative dx/dt(t) at a certain time t I need the results of x from t0 to t (t0<t). So the computation requires the former results of x. Someone suggested to me, that I can use Simulink with a "Level 2 S function"-block. Within the Setup-function of the "Level 2 function" I use the "Dwork" element to store the time history of the required parts of x (I only need the history of x(1) and x(9) ). The results shall be stored every second (the dynamic is quite slow). FzSeedStateModel returns the derivative. The history stored in "Dwork" is forwarded by the local function

param.seedShape=@(s_s) getSeedShape(block,s_s);

My currently (not working) code:

function dynamic2_l2fct(block)
setup(block);

%endfunction

%%Function: setup ===================================================
%%Abstract:

Set up the basic characteristics of the S-function block such as:

- Input ports

- Output ports

- Dialog parameters

- Options

Required : Yes

C-Mex counterpart: mdlInitializeSizes

function setup(block)

% Register number of ports
block.NumInputPorts  = 1;
block.NumOutputPorts = 1;

% Setup port properties to be inherited or dynamic
block.SetPreCompInpPortInfoToDynamic;
block.SetPreCompOutPortInfoToDynamic;

% Override input port properties
% input:
%   v_s         seed pull velocity
%   v_f         feed rod push velocity
%   P           heating power
block.InputPort(1).Dimensions        = 3;
block.InputPort(1).DatatypeID  = 0;  % double
block.InputPort(1).Complexity  = 'Real';
block.InputPort(1).DirectFeedthrough = false;

% Override output port properties
block.OutputPort(1).Dimensions       = 10;
block.OutputPort(1).DatatypeID  = 0; % double
block.OutputPort(1).Complexity  = 'Real';

% Register parameters
% parameters are:
%   param       containing physical parameters
%   x0          initial state
block.NumDialogPrms     = 2;

% Set up the continuous states.
block.NumContStates = 1;

% Register sample times
%  [0 offset]            : Continuous sample time
%  [positive_num offset] : Discrete sample time
%
%  [-1, 0]               : Inherited sample time
%  [-2, 0]               : Variable sample time
block.SampleTimes = [0 0];  % for continuous states
% block.SampleTimes = [1 0];    % for discrete states

% Specify the block simStateCompliance. The allowed values are:
%    'UnknownSimState', < The default setting; warn and assume DefaultSimState
%    'DefaultSimState', < Same sim state as a built-in block
%    'HasNoSimState',   < No sim state
%    'CustomSimState',  < Has GetSimState and SetSimState methods
%    'DisallowSimState' < Error out when saving or restoring the model sim state
block.SimStateCompliance = 'DefaultSimState';

%%-----------------------------------------------------------------
%%The MATLAB S-function uses an internal registry for all
%%block methods. You should register all relevant methods
%%(optional and required) as illustrated below. You may choose
%%any suitable name for the methods and implement these methods
%%as local functions within the same file. See comments
%%provided for each function for more information.
%%-----------------------------------------------------------------

block.RegBlockMethod('PostPropagationSetup',    @DoPostPropSetup);
block.RegBlockMethod('InitializeConditions', @InitializeConditions);
block.RegBlockMethod('Start', @Start);
block.RegBlockMethod('Outputs', @Outputs);     % Required
block.RegBlockMethod('Update', @Update);        % no discrete states
block.RegBlockMethod('Derivatives', @Derivatives);
% block.RegBlockMethod('Terminate', @Terminate); % Required

%end setup
%%PostPropagationSetup:

Functionality : Setup work areas and state variables. Can

also register run-time methods here

Required : No

C-Mex counterpart: mdlSetWorkWidths

function DoPostPropSetup(block)
block.NumDworks = 3;
% block.NumContStates = 1;
  % state for seed length
  block.Dwork(1).Name            = 'l_s';
  block.Dwork(1).Dimensions      = 1e4;
  block.Dwork(1).DatatypeID      = 0;      % double
  block.Dwork(1).Complexity      = 'Real'; % real
  block.Dwork(1).UsedAsDiscState = true;
  % state for seed radius
  block.Dwork(2).Name            = 'R_s';
  block.Dwork(2).Dimensions      = 1e4;
  block.Dwork(2).DatatypeID      = 0;      % double
  block.Dwork(2).Complexity      = 'Real'; % real
  block.Dwork(2).UsedAsDiscState = true;
  % state to store the current index (number loop iterations)
  block.Dwork(3).Name            = 'idx';
  block.Dwork(3).Dimensions      = 1;
  block.Dwork(3).DatatypeID      = 0;      % double
  block.Dwork(3).Complexity      = 'Real'; % real
  block.Dwork(3).UsedAsDiscState = true;
  % continuous state
%   block.ContStates.Name            = 'state';
%   block.ContStates.Dimensions      = 10;
%   block.ContStates.DatatypeID      = 0;      % double
%   block.ContStates.Complexity      = 'Real'; % real
%   block.ContStates.UsedAsDiscState = false;
%%InitializeConditions:

Functionality : Called at the start of simulation and if it is

present in an enabled subsystem configured to reset

states, it will be called when the enabled subsystem

restarts execution to reset the states.

Required : No

C-MEX counterpart: mdlInitializeConditions

function InitializeConditions(block)
    % assign initial states
    block.ContStates.Data=block.DialogPrm(2).Data;
    % initial seed length (usually zero)
    block.Dwork(1).Data(1)=block.DialogPrm(2).Data(9);
    % initial seed radius
    block.Dwork(2).Data(1)=block.DialogPrm(2).Data(1);
    % first entry in Dwork(i).Data(1) is already done
    % Data(2) will be next
    block.Dwork(3).Data=2;
%end InitializeConditions
%%Start:

Functionality : Called once at start of model execution. If you

have states that should be initialized once, this

is the place to do it.

Required : No

C-MEX counterpart: mdlStart

function Start(block)

% block.Dwork(1).Data = 0;
%end Start
%%Outputs:

Functionality : Called to generate block outputs in

simulation step

Required : Yes

C-MEX counterpart: mdlOutputs

function Outputs(block)

%   block.OutputPort(1).Data = block.Dwork(1).Data + block.InputPort(1).Data;
    block.OutputPort(1).Data = block.ContStates(1).Data;

%end Outputs
%%Update:

Functionality : Called to update discrete states

during simulation step

Required : No

C-MEX counterpart: mdlUpdate

function Update(block)
i=block.Dwork(3).Data;
block.Dwork(1).Data(i) = block.ContStates(1).Data(9);
block.Dwork(2).Data(i) = block.ContStates(1).Data(1);
block.Dwork(3).Data = i+1;
%end Update
%%Derivatives:

Functionality : Called to update derivatives of

continuous states during simulation step

Required : No

C-MEX counterpart: mdlDerivatives

function Derivatives(block)
    % need to create function and fuction handle to compute the seed shape
    param=block.DialogPrm(1).Data;
    param.seedShape=@(s_s) getSeedShape(block,s_s);
    block.Derivatives(1).Data=FzSeedStateModel(block.CurrentTime,block.ContStates(1).Data,block.InputPort(1).Data,param);
%end Derivatives
%%Terminate:

Functionality : Called at the end of simulation for cleanup

Required : Yes

C-MEX counterpart: mdlTerminate

function Terminate(block)

%end Terminate

function varargout=getSeedShape(block,s)
    % check, if s is in (-1,0) 
    if(any(s<-1) || any(s>0))
        error('Error in seed arc length!')
    end
    % convert arc length to length argument (l_s+z, + because z<0)
    z=block.DialogPrm(1).Data.Z_s.*s; % s in [-1,0] -> z in [-Z_s,0]
    l=block.ContStates(1).Data(9)+z; % z in [-Z_s,0] -> l in [l_s, l_s-Z_s]
    % now find the appropriate radius
    r=zeros(size(l));
    for j=1:length(l)
        % if l(j) < 0 -> seed isn't long enough, take radius obtained from
        % steady state seed shape function
        if(l(j)<0)
            r(j)=block.DialogPrm(1).Data.seedShape(l(j),block.DialogPrm(1).Data,'height');
        end
        % index of last length l_s that is lager than the length l(j)
        % last index because this is the "youngest"
        idx=find(block.Dwork(1).Data>l(j),1,'last');
        if(isempty(idx))
            % seed isn't long enouth
            % -> take initial seed radius
            r(j)=block.DialogPrm(2).Data(1);
        else
            % seed is long enough
            % check if previous length is less than l(j) otherwise error
            if(l(j)<block.Dwork(1).Data(idx-1))
                error('error when getting correct length index');
            end
            % interpolate data to get radius
            r(j)=interp1(block.Dwork(1).Data(idx-1:idx),block.Dwork(2).Data(idx-1:idx),l(j));
        end
    end

      if(nargout==0 || nargout==1)
          varargout{1}=[r,z];
      elseif(nargout==2)
          varargout{1}=r;
          varargout{2}=z;
      else
          error('Output for number of required outputs not defined!');
      end

Currently, I get the error that I want to assign a vector of size ten to a vector of size 1. But I don't know where to define the dimension of the continuous state. In the "Setup" function it s prohibited and later (in "Start" or "InitializeConditions") it fails.

My questions are: -

• is it possible to define one continuous state with an dimension of ten? (something like block.ContStates.Dimension=10)? And if this is possible, where and how to define this? Or does 10 ODE's require 10 continuous states?
• how to achieve, that "Update" will be executed only every second, whereas the ODE's are integrated by variable step solver (i.e. ode45)

Or is there an other possibility to integrate a system of ODE's using the former results?

Please excuse my bad english.