I had skipped an important step in my implementation and a few other things were wrong too. After fixing those issues and playing around with the parameters a bit i found a working combination so i can give you an example.
Let p = 289254655033, q = 576650390627 and N = p*q = 166798809815457581575691. We want to factor N. We setup some parameter which determine the polynomial f. By construction the polynomial f has a root (x_0, y_0) which will give us the factorization of N. It holds that p = m*x_0 + c and q = n*y_0 + d where we have to determine c,d by brute force search and x_0,y_0 by finding an appropriate second polynomial g with the same root. In this example is m = 5779 and n = 4877 (determined by my parameter choice). We want both m and n to be as small as possible because their size determines how many combinations for c and d we have to test. By brute force search (running the entire algorithm again and again) we find c = 3490 and d = 2000.
We obtain the polynomial f( x, y ) = 272663114441221050502036905165359*x*y + 111815917343129403527593563734000*x + 164664175359034688744092195713290*y - 1613666891439484824286517533623583653166443925443
Now we construct a lattice and use lattice reduction techniques to find a short coefficient vector for our new polynomial g. This polynomial must have the same root and may not be a multiple of f (gcd(f,g) = 1 because f is irreducible).
We find: g( x,y ) = -550027317057214807707838758454131*x^2*y^2 - 21035627259174470964704774394001486176978740190*x^2*y + 279798978942801467799047317444099858865221454226*x*y^2 - 8626461865562536561414577174664878412692940000*x^2 + 815210635045419056587043555126594181612580850722711584*x*y + 168973600365180529829184985202345054933202887160*y^2 - 138762500323450558676026989874427343055634473948375906415055203127687897949926*x + 60213348692682732919411494951824394082422924853210355495670640666282857023814*y - 174110984048404287435752616687571392254879687559395285145418320554729063981835001233
Now we compute the resultant of f and g:
4476571884282525266836510330141869220498240713975097720401448926613563671146029895483060633149767902954109339351810218406344618956582075732396*y^3 - 12944301498631912946627224525490850239753133562803345399593047007314825269294562029498487898880328775532821414009302732104707106625214703423125443824*y^2 - 61053757336251582441768270515159399095157293517551804381202990288793114577733640928959281804253883501665688671801616653381775535847520521818182191534537076572*y - 25037423551750369812975300235336096113581272995509465458164313917743589489401221657131723392563075939670940390415678987548009123529431653353682394713202072000
The roots of this polynomial are y_0 = 118238751 and y_0 = 1. We can disregard y_0 = 1 (don't know where this comes from) and use the special construction of f to obtain q = 4877*118238751+2000 = 576650390627. So we are done.
The implementation:
in sage
Note: I still don't use brute force search here and instead calculate c = p mod m and d = q mod n so don't be confused. Also the calculation takes place over the integers obviously.