Like scipy.optimize.fsolve, scipy.optimize.root finds the solution of a nonlinear vector valued function, that is, the the solution of a system of nonlinear equations. Based on this it appears that the difference between root and fsolve is that fsolve requires the dimension of the input (e.g. x0) to be the same as the dimension of the output (e.g. return of f) but I haven't confirmed by testing. I agree that this is not immediately apparent from the documentation (and this might be something you could clear up). Also, can others chime in on why they are separate? Backwards compatibility? 

Matt