You use the regress() method for each least squares fit required by the algoritm. In essence, given three lines of regression from each plane in the model, you treat those as equations and use a matrix approach to solving them as simultaneous equations to produce the resulting plane (PM supersearch for "zentara" and "solving matrices" for a more detailed explanation of that sub-algorithm) - there isn't a method that does the whole algorithm for you.