function [g]=gsystem(t,y,p,xdim,kdim,solution_only)
%
%   ***   function [g]=gsystem(t,y,p,xdim,kdim,solution_only)   ***
%
%  Input
%     t = time
%     y = current solution of length xdim*(1+kdim+xdim) (=ny)
%     p = vector of parameters of length kdim
%     xdim = number of variables
%     kdim = number of parameters
%     solution_only = Compute only the solutions and not derivatives
%
%  Output
%      g = rhs vector for the system of odes of length xdim*(1+kdim+xdim)

% Establish g and dgdp
ny=length(y);
g=zeros(ny,1);
dgdp=zeros(xdim,kdim+xdim);

% Establish x and dxdp
x=zeros(xdim,1);
dxdp=zeros(xdim,kdim+xdim);

% Unpack x vector from y vector
x=y(1:xdim);

% Unpack dxdp matrix from y vector
if solution_only == 0
    dxdp=reshape(y(xdim+1:ny), xdim, kdim+xdim);
end

% Calculate g
g(1:xdim)=gvec(x,p);

% Calculate dgdp
if solution_only == 0
    dgdp=drhs_dparam(x,dxdp,p);

    % Repack dgdp as a column vector and add to g
    % (In fact this is unnecessary as it is the default, but it makes the point)
    g(xdim+1:ny) =reshape(dgdp,xdim*(kdim+xdim),1);
end

end