PHASEGO: a toolkit for automatic calculation and plot of phase diagram. The PHASEGO package extracts the Helmholtz free energy from the phonon density of states obtained by the first-principles calculations. With the help of equation of states fitting, it reduces the Gibbs free energy as a function of pressure/temperature at fixed temperature/pressure. Based on the quasi-harmonic approximation (QHA), it calculates the possible phase boundaries among all the structures of interest and finally plots the phase diagram automatically. For the single phase analysis, PHASEGO can numerically derive many properties, such as the thermal expansion coefficients, the bulk moduli, the heat capacities, the thermal pressures, the Hugoniot pressure-volume-temperature relations, the Gruneisen parameters, and the Debye temperatures. In order to check its ability of phase transition analysis, I present here two examples: semiconductor GaN and metallic Fe. In the case of GaN, PHASEGO automatically determined and plotted the phase boundaries among the provided zinc blende (ZB), wurtzite (WZ) and rocksalt (RS) structures. In the case of Fe, the results indicate that at high temperature the electronic thermal excitation free energy corrections considerably alter the phase boundaries among the body-centered cubic (bcc), face-centered cubic (fcc) and hexagonal close-packed (hcp) structures.