The paper puts forward an effective algorithm for producing approximate polynomial solutions for linear ordinary differential equations (LODEs) and sets of LODEs with polynomial coefficients and polynomial right-hand side functions. The algorithm is an upgraded version of the Lanczos Tau-method and provides the optimal deviation of the approximate solution from the exact one according to the minimax norm for a given interval. With minor modification, the algorithm allows one to find approximate expressions for the derivatives of the exact solutions with sufficiently greater accuracy than the derivatives of the approximate solutions are capable of providing that.