function [dxdt Jx Jp] = kron2eqns(m)

eqns = cell(m.nX, 1);

% write equations out for the a1P matrix
[i j] = find(m.a1P);

ix = mod(i-1, m.nX) + 1;
ip = floor((i-1)/m.nX) + 1;

for r = 1:length(ix)
    if isempty(eqns{ix(r)})
        eqns{ix(r)} = ['ddt_' m.xNames{ix(r)} ' = ' eqns{ix(r)} num2str(m.a1P(i(r), j(r))) '*' m.pNames{ip(r)} '*' m.xNames{j(r)}];
    else
        eqns{ix(r)} = [eqns{ix(r)} ' + ' num2str(m.a1P(i(r), j(r))) '*' m.pNames{ip(r)} '*' m.xNames{j(r)}];
    end
end

% write equations out for the a2P matrix

[i j] = find(m.a2P);

ix = mod(i-1, m.nX) + 1;
ip = floor((i-1)/m.nX) + 1;
ix1 = mod(j-1, m.nX) + 1;
ix2 = floor((j-1)/m.nX) + 1;

for r = 1:length(ix)
    if isempty(eqns{ix(r)})
        eqns{ix(r)} = ['ddt_' m.xNames{ix(r)} ' = ' eqns{ix(r)} num2str(m.a2P(i(r), j(r))) '*' m.pNames{ip(r)} '*' m.xNames{ix1(r)} '*' m.xNames{ix2(r)}];
    else
        eqns{ix(r)} = [eqns{ix(r)} ' + ' num2str(m.a2P(i(r), j(r))) '*' m.pNames{ip(r)} '*' m.xNames{ix1(r)} '*' m.xNames{ix2(r)}];
    end
end

% write equations out for the b1P matrix ** NOT YET DEBUGGED! **

[i j] = find(m.b1P);

ix = mod(i-1, m.nX) + 1;
ip = floor((i-1)/m.nX) + 1;

for r = 1:length(ix)
    if isempty(eqns{ix(r)})
        eqns{ix(r)} = ['ddt_' m.xNames{ix(r)} ' = ' eqns{ix(r)} num2str(m.b1P(i(r), j(r))) '*' m.pNames{ip(r)} '*' m.uNames{j(r)}];
    else
        eqns{ix(r)} = [eqns{ix(r)} ' + ' num2str(m.b1P(i(r), j(r))) '*' m.pNames{ip(r)} '*' m.uNames{j(r)}];
    end
end

% write equations out for the b2P matrix
[i j] = find(m.b2P);

ix = mod(i-1, m.nX) + 1;
ip = floor((i-1)/m.nX) + 1;
ix1 = mod(j-1, m.nX) + 1;
iu1 = floor((j-1)/m.nX) + 1;

for r = 1:length(ix)
    if isempty(eqns{ix(r)})
        eqns{ix(r)} = ['ddt_' m.xNames{ix(r)} ' = ' eqns{ix(r)} num2str(m.b2P(i(r), j(r))) '*' m.pNames{ip(r)} '*' m.xNames{ix1(r)} '*' m.uNames{iu1(r)}];
    else
        eqns{ix(r)} = [eqns{ix(r)} ' + ' num2str(m.b2P(i(r), j(r))) '*' m.pNames{ip(r)} '*' m.xNames{ix1(r)} '*' m.uNames{iu1(r)}];
    end
end


% write equations out for the kP vector ** NOT YET DEBUGGED **

i = find(m.kP);
ix = mod(i-1, m.nX) + 1;
ip = floor((i-1)/m.nX) + 1;

for r = 1:length(ix)
    if isempty(eqns{ix(r)})
        eqns{ix(r)} = ['ddt_' m.xNames{ix(r)} ' = ' eqns{ix(r)} num2str(m.kP(i(r))) '*' m.pNames{ip(r)}];
    else
        eqns{ix(r)} = [eqns{ix(r)} ' + ' num2str(m.kP(i(r))) '*' m.pNames{ip(r)}];
    end
end


% this line replaces complexes delimited by the ':' character with the '_'
% character. because this is commmon parlance in Kronecker and the
% character cannot be used in matlab variables, this makes the translation
% process easier!
eqns = strrep(eqns, ':', '_');

% % print output if necessary
% for i = 1:m.nX
%     fprintf(1, '%s\n', eqns{i});
% end

%% symbolic manipulation (simplification) and Jacobian computation

for i = 1:m.nX
    syms(strrep(m.xNames{i}, ':', '_'));
end

for i = 1:m.nP
    syms(m.pNames{i});
end

for i = 1:m.nU
    syms(m.uNames{i});
end

dxdt = cell(m.nX, 1);
Jx   = cell(m.nX, m.nX);
Jp   = cell(m.nX, m.nP);

for i = 1:length(eqns)
    eval([eqns{i} ';']);
    dxdt{i} = char(eval(['ddt_' strrep(m.xNames{i}, ':', '_')]));
    for j = 1:m.nX
        Jx{i,j} = char(eval(['diff(' ['ddt_' strrep(m.xNames{i}, ':', '_') ] ', ' strrep(m.xNames{j}, ':', '_') ');']));
    end
    for j = 1:m.nP
        Jp{i,j} = char(eval(['diff(' ['ddt_' strrep(m.xNames{i}, ':', '_') ] ', ' strrep(m.pNames{j}, ':', '_') ');']));
    end
end