diff --git a/src/quadratization.jl b/src/quadratization.jl
new file mode 100644
index 0000000..acf4ea7
--- /dev/null
+++ b/src/quadratization.jl
@@ -0,0 +1,49 @@
+"""
+Abstract type for quadratization algorithms.
+"""
+abstract type QuadratizationAlgorithm end
+
+"""
+    QBee <: QuadratizationAlgorithm
+
+QBee (Quadratization via monomial introduction) algorithm from Bychkov & Pogudin (2021).
+
+This algorithm finds an optimal monomial quadratization by introducing new variables
+that represent monomials in the original state variables. It minimizes the number of
+new variables needed to transform the system to quadratic form.
+
+# Reference
+Bychkov, A., & Pogudin, G. (2021). Optimal monomial quadratization for ODE systems.
+https://arxiv.org/abs/2103.08013
+"""
+struct QBee <: QuadratizationAlgorithm end
+
+"""
+    quadratize(sys::System, alg::QuadratizationAlgorithm = QBee(); 
+               name = Symbol(nameof(sys), :_quadratized), kwargs...)
+
+Transform a polynomial ODE system into quadratic form using the specified algorithm.
+
+# Arguments
+- `sys::System`: Input polynomial ODE system with first-order derivatives on LHS
+- `alg::QuadratizationAlgorithm`: Algorithm to use for quadratization (default: `QBee()`)
+- `name`: Name for the output system (default: name of input system)
+- `kwargs...`: Algorithm-specific keyword arguments
+
+# Returns
+- `System`: Quadratic ODE system with introduced auxiliary variables
+"""
+function quadratize(
+        sys::System, alg::QuadratizationAlgorithm = QBee();
+        name = Symbol(nameof(sys), :_quadratized), kwargs...
+    )::System
+    return quadratize(sys, alg, name; kwargs...)
+end
+function quadratize(sys::System, alg::QBee, name; kwargs...)::System
+    return quadratize_qbee(sys, name; kwargs...)
+end
+
+function quadratize_qbee(sys::System, name; kwargs...)
+    # TODO
+    return sys
+end
diff --git a/src/utils.jl b/src/utils.jl
index d134175..a10c73d 100644
--- a/src/utils.jl
+++ b/src/utils.jl
@@ -126,3 +126,26 @@ function separate_terms(exprs::AbstractVector, vars, iv)
     linear = sparse(linear_I, linear_J, linear_V, length(exprs), length(vars))
     return linear, other_terms, nonlinear_terms
 end
+
+"""
+$(TYPEDSIGNATURES)
+
+Convert a symbolic expression to an expanded polynomial representation of
+`DynamicPolynomials.Polynomial` monomials.
+
+This function transforms symbolic polynomial expressions into a monomial form that can be 
+more easily processed in certain applications.
+"""
+function to_expanded_polynomial(expr::Num)
+    val = Symbolics.unwrap(expr)
+    return to_expanded_polynomial(val)
+end
+function to_expanded_polynomial(expr::SymbolicUtils.BasicSymbolic{T})::
+    Union{SymbolicUtils.PolyVarT, SymbolicUtils.PolynomialT} where {T}
+    if !SymbolicUtils.iscall(expr)
+        return expr
+    end
+    poly_to_bs = Dict{SymbolicUtils.PolyVarT, SymbolicUtils.BasicSymbolic{T}}()
+    bs_to_poly = Dict{SymbolicUtils.BasicSymbolic{T}, SymbolicUtils.PolyVarT}()
+    return SymbolicUtils.to_poly!(poly_to_bs, bs_to_poly, expr)
+end