%
 function M_dot = rhs_post(t,M,FLAG,local,data)
 %
 % M is a concatention of  two matrices
 % M = [\Phi, \Gamma], where \Phi is a solution of the matrix equation
 % \Phi_dot = J_x \Phi
 % with initial condition \Phi(t_l ) = eye.
 % J_x  is the (state) Jacobian \frac{\partial f]{\partial x}.
 % \Gamma is the solution of the matrix equation
 % \Gamma_dot = J_x  \Gamma +  J_u  
 % with initial condition \Gamma(t_l) = zeros(N_states,N_controls).
 %  J_u  is the (control) Jacobian \frac{\partial f]{\partial u}.
 %
     [ J_x, J_u ] = jacobian(t,local,data.Rhsj);
 %
     Phi_dot = J_x*reshape(M(1:data.N_states*data.N_states), ...
                                     data.N_states,data.N_states);
     Gam_dot = J_x*reshape(M(1+data.N_states*data.N_states:end), ...
                          data.N_states,data.N_controls) + J_u;
%
    M_dot = [Phi_dot Gam_dot];  
    M_dot = M_dot(:);
%
% return
%
%========================================================== 
% internal functions for rhs_post
%========================================================== 
%
 function [ J_x, J_u ] = jacobian(t,loc,rhsj)
%
%  J_x  is the (state)   Jacobian \frac{\partial f]{\partial x}.
%  J_u  is the (control) Jacobian \frac{\partial f]{\partial u}.
% 
% local is a data structure that has 
%
%   loc.t_l  the time value at the 'left' end of this interval
%   loc.t_r  the time value at the 'right' end of this interval
%   loc.x_l  the state value at the 'left' end of this interval
%   loc.t_r  the state value at the 'right' end of this interval
%   loc.u    the control vector on this interval
%
%   rhsj      the name of the routine for calculating the Jacobian matrices
%
% first interpolate to estimate x at this time
%
    if ( t > loc.t_r  |  t < loc.t_l)
        t
        loc.t_l
        loc.t_r
        error(['error in computing the Jacobian'])
    end
%
   x = ((loc.t_r - t)*loc.x_l + (t - loc.t_l)* loc.x_r)/(loc.t_r -loc.t_l);
%
   [ J_x, J_u ]  = feval(rhsj,t,x,loc.u);
%
% return
%

