We present a solution of the DGLAP evolution equations,
written in terms of Sudakov form factors to describe the branching
and no-branching probabilities, using a Monte Carlo parton branching method. We
demonstrate numerically that this method reproduces
the semi-analytical solutions.
We show how this method can be used to determine
Transverse Momentum Dependent (TMD) parton distribution functions,
in addition to the usual integrated parton distributions functions.
We discuss numerical effects of the boundary of soft gluon resolution scale parameter on the resulting parton distribution functions.
We show that a very good fit of the integrated TMDs to high precision HERA data can be
obtained over a large range in $x$ and $Q^2$.