Superintegrable and shape invariant systems with position dependent mass

Second order integrals of motion for 3d quantum mechanical systems with position dependent masses (PDM) are classified. Namely, all PDM systems are specified which, in addition to their rotation invariance, admit at least 1 second order integral of motion. All such systems appear to be also shape invariant and exactly solvable. Moreover, some of them possess the property of double shape invariance and can be solved using two different superpotentials. Among them there are systems with double shape invariance which present nice bridges between the Coulomb and isotropic oscillator systems. A simple algorithm for calculating the discrete spectrum and the corresponding state vectors for the considered PDM systems is presented and applied to solve five of the found systems.