circuit simulation involving nonlinear elements can be a challenging task. from these challenges rises a demand for finding better and more efficient ways of solving such problems. this work provides a novel approach for performing harmonic balance (hb) analysis using the freeda circuit simulator. the proposed method is an extension of the method of multiple reflections for multiple ports. in addition the method is formulated in terms of power waves and state variables at the nonlinear devices. the hb problem is then solved using a procedure which resembles the signal propagation within the actual circuit. this method could be efficiently parallelized since it does not require a large matrix decompositions at each iteration. several approaches to improve convergence properties are investigated. the first involves adding capacitors in parallel with the nonlinear device ports, this allows the fixed-point iterations to always be convergent. these capacitors are only active in a separate time dimension and do not affect the steady-state solution. the harmonic balance solution is found when the transient response in this time dimension is extinguished. another strategy to improve convergence is the combination of fixed-point iterations with the gradient descent method. the effect of a vector extrapolation method to accelerate convergence is also investigated. simulation results for various strongly nonlinear circuits is presented. this thesis covers the background of harmonic balance analysis, literature review, derivation of the proposed method, improvements, preliminary results, as well as future work.