You will need to write a solver for the incomplete combustion of your propellant. This can be done by hand, but I recommend writing a program to implement it.
These days, such calculations are primarily done by solving for the maximum entropic state of the system (the system being the mixture of reactants at the chamber temperature & pressure). The equivalent problem to the maximization of the entropic state is the minimization of the Gibbs free energy (or the Helmholtz free energy, depending on whether you allow pressure or temperature, respectively, to vary in your calculations).
Traditionally, Gibbs minimization problems are solved in R-space, where R is the number of different reactions between species. Because we are modeling non-ideal combustion, we do not just consider the dominant species. For example, consider the very simple reaction 2H2 + O2, a hydrolox reaction.
Besides the obvious H2O product, you must also consider H2, O2, H, O, O3, OH, etc. and, depending on the robustness of your method, e- and all associated molecular and monatomic ions. Multiple reactions might be considered per species. For example, 2H2 + O2 --> 2H2O, but so does OH + H --> H2O, O + H2 --> H2O, OH- + H+ --> H2O, H3O+ + e- --> H2O + H, etc. All reverse equations are also implicitly considered during solving. As you can see, this produces a very large, non-linear system. The hydrolox reaction is the simplest you are likely to find in rocketry. Large CHON systems can grow to hundreds of reactions, & thus dimensions.
However, by the method of Reynolds (1986), you may reduce the dimensionality by rejecting the traditional method of balancing the reactions & instead converting the problem to a less-physical one. The concept of element potentials arise when implementing the constraint of atomic population (aka: no matter how the atoms are redistributed in the product species, there must be no net change in the total number of atoms in the system, per change in Gibbs energy) as Lagrangian multipliers.
As these Lagrange multipliers are as many as the types of atom in the system, they far less numerous than the number of reactions, or even the number of species. A CHON system only has 4 types of element, but it can have 100s of reactions & dozens of species. Remarkably, the element potentials can be solved for, and in doing so, solve the system. From the aforementioned paper:
This method is also more robust to the influence of rare but high-energy species, such as solid phases like elemental C (soot), and ionized species. These species can have a large influence on the thermodynamics of the system while existing in population several orders of magnitude lower than the reactants & dominant products.
While you are still solving a non-linear system, Reynolds demonstrates that the dual problem can be solved by a minimaxing method of steepest descent/ascent, followed by the applications of a second-order Newton method to rapidly approach the optimal solution from near the solution, where the solution space shallows out, becomes difficult for gradient descent methods to operate, and the space becomes well-approximated by second-order polynomials.
I am still parsing the work of Gordon & McBride (1994) but I will update my answer with their additions to the science when I am done re-reading their paper for the umpteenth time.
Once you have the energy of your reaction, it is relatively trivial to determine your thrust & specific energy using the adiabatic & isentropic mass flow equations. You will also know the composition of your exhaust, so you can determine the parameter gamma, the ratio of specific heats, which is required for said mass flow equations.
