It's quite important that you create a model as accurately as possible, both for making your artifact and debugging your code. For example, if you specify an axis-aligned rectangular polygon, but the 4 corners don't lie exactly on the axes, when you compute the 4 corners in turn they may end up not lying on the same plane.
Instead of simply judging with your eyes and placing points "roughly" where they should go, use the Guidelines feature so the program can help you place points with perfect accuracy. See the UI page for a description on how to use this feature.
Another tip is to first specify points as described above, then use them as the corners of your polygon, instead of directly creating a polygon from scratch. Points can be reused for new lines or polygons by holding down the Ctrl key and clicking on them when you specify lines or polygons.
A third tip is to work on a closely zoomed in image, so you can get subpixel accuracy. You can zoom in and out on an image by Ctrl-+/- (no need to hold down Shift for +).
Note that you can still get a small amount of error due to certain inherent limitations, such as the vanishing points (they are the "best" fit, not an exact fit), the need to specify points where Guideline feature cannot help, etc.
Although this document describes how to create the plane coordinate system, it doesn't explain how to choose the p, q and r.
Unfortunately you can't simply choose the first 3 corners of the polygon as your p, q, r, because you don't want all 3 of them to lie on the same line. And this is quite possible because you may have added an extra collinear point by mistake (e.g. clicking on the first corner again before pressing enter to finish the polygon, which effectively adds another point very close to the first point) or intentionally (to use as the known point for SameXY or SameZPlane calculation).
One way to ensure that you choose appropriate p, q, r is following (which the sample solution does). Set the 1st and 2nd corner as your p and r, respectively. Then search among the rest of points the point X for which the angle between (p-r) and (X-r) is closest to 90 degrees, and use that as your q.
This can be achieved by looking at the dot product of normalized (p-r) and normalized (X-r), which is equivalent to the cosine of the angle between them, and choose X whose cosine is closest 0.
As described in this write-up, the line endpoints are translated (u and v) and scaled (w) before passing them into the matrix. Once you obtain the eigenvector, you must make sure to apply the inverse of the above transformations so that the eigenvector has the correct image coordinate. And it's important that they are applied in the right order.
First, scale the *entire* eigenvector so that its 3rd element (corresponding to w) is the same as when it was scaled before matrix computation. Then apply the inverse translation to the 1st and 2nd element (corresponding to u and v) to obtain the final intersection point.