function [TdotXdot] = TdotXdot_func( t, TX, L, q_M )
% function to compute the derivative of temperature at all
%  the nodal locations in a 1-D spatially discretized finite 
%  element domain.  We will use a Runge-Kutta like 
%  scheme to perform the time integration.
%  The derivatives of the sensitivity problem are also computed
%  and appended to the Tdot vector (T and X must be computed
%  in parallel due to the dependence on the thermal properties
%
%if( nargin < 4 ) 
%    nargin
%    error('Must provide L and q_M in Tdot');
%end
%
n_nodes = length( TX )/2;  % how many are there
T = TX(1:n_nodes);         % get T only
X = TX(n_nodes+1 : 2*n_nodes );  % this is X only
n_ele = n_nodes - 1;    % for 1-D there are one less elements
dx = L / n_ele;
for i = 1:n_ele
    Tavg(i,1) = (T(i)+T(i+1))/2;
end
Cvec = C_func( Tavg );  % compute the element values
kvec = k_func ( Tavg );

% assemble the conductivity and capacitance matricies
K_c = zeros( n_nodes, n_nodes );
C = zeros( n_nodes, n_nodes );
for i = 1:n_ele
    K_c(i,i) = K_c(i,i) + kvec(i)/dx;
    K_c(i,i+1) = K_c(i,i+1) - kvec(i)/dx;
    K_c(i+1,i) = K_c(i+1,i) - kvec(i)/dx;
    K_c(i+1,i+1) = K_c(i+1,i+1) + kvec(i)/dx;
    C(i,i) = C(i,i)+ Cvec(i)*dx/3;
    C(i+1,i+1) = C(i+1,i+1)+ Cvec(i)*dx/3;
    C(i+1,i) = C(i+1,i)+ Cvec(i)*dx/6;
    C(i,i+1) = C(i,i+1)+ Cvec(i)*dx/6 ;
end
% compute the RHS vector - only unknown is q_M
RHS = zeros( n_nodes, 1);
RHS(1) = q_M;
% RHS(1) = q_func( t );
Tdot = C \ ( RHS - K_c * T );
RHS(1) = 1.0;    % for the Xdot problem
Xdot = C \ ( RHS - K_c * X );
TdotXdot = [ Tdot; Xdot ];

return

function [C] = C_func( T )
% C = ones( length( T ) );
% C = 2.5 ./ ( T + 0.5 );
C = 1 ./ ( 2*T + 1 );
return

function [k]=k_func(T)
% k = ones( length(T) );
% k = -25*T.*T + 20 * T + 1;
k = 1+4*T.^2;
return


function [q] = q_func( t )
if( t < 0 )
    q = 0;
elseif( t < 0.6 )
    q = t;
elseif( t < 1.2 )
    q = 1.2 - t;
else
    q = 0;
end
return
