ddt(of, n,state,dstate,f,u)

function DNof = ddt(of, n,state,dstate,f,u)
%DDT       Calculate time derivative as a function of state and input variables of a (nonlinear) dynamic system
% function ddt(of, n, state,dstate,f,u)
%
% Given:
%   of:     the symbolic function you want the derivative of (e.g. x2)
%   n:      n'th order derivative
%   state:  a column vector denoting the names of the state, e.g. [x1;x2]
%   dstate: the names of the states plus a d in front of them, e.g. [dx1; dx2]
%   f:      symbolic expression for the derivative of the state, e.g. [x2, sin(u1) ]
%   u:      a list of inputs, e.g. [u1]
%
% Example:
%   Suppose we have the following (nonlinear) dynamic system:
%
%       d/  / x1 \    / 2 * cos(x2) + u1 \
%       dt  \ x2 / =  \  - x1^2 + x1*u2  /   , y = x2 + x1 + u2
%
%   where x1, x2 are time-dependent states and u1/u2 are time-dependent inputs
%   and we are interested (for some reason) in the second time derivative of
%   y, expressed in terms of the state [x1;x2] and the inputs and its
%   derivatives ([u1;u2], [du1;du2], [d2u1; d2u2]). We calculate this as
%   follows:
%
%     clear
%     syms x1 x2 dx1 dx2 u1 u2;
%     eqn = [ 2*cos(x2) + u1; -x1^2 + x1*u2];
%     y = x2 + x1 + u2;
%     ddy = ddt ( y, 2, [x1;x2], [dx1;dx2], eqn, [u1;u2])
%
%   Output:
%     ddy =
%        (-2*x1+u2)*(2*cos(x2)+u1)-2*sin(x2)*(-x1^2+x1*u2)+du1+x1*du2+d2u2
%
% Note: This function needs the function 'fulldiff' by Tim Jorris, which
%   can be downloaded from Matlab's File Exchange.

% History:
% 081016 (GvO) Added example in help, Put into SVN
% 07xxxx (GvO) Made
for (i = 1:n) % do n derivations
    timedependents = num2cell([state;u]);
    dof = fulldiff(of,timedependents);
    of = subs(dof,dstate,f);
    sof = simple(of);
end
DNof = sof;