The present study investigates the thermoelastic response of a heterogeneous piezoelectric sphere under thermal shock loading. Boundary conditions as well as loading are considered as symmetric; thus, the response of the sphere is expected to be symmetric. All of the properties of the thick-walled sphere, including mechanical, electrical and thermal properties, are considered dependent on the radial position, except for the relaxation time, which is considered a constant value along the radius. The governing equations of the sphere have been derived under heterogeneous anisotropic assumptions. The general form of the second law of thermodynamics, which is nonlinear in nature, and is called nonlinear energy equation is used. The number of the established equations is three, which includes the motion equation, energy equation and Maxwell electrostatic equation of Maxwell. These equations are obtained in terms of radial displacement, temperature difference and electric potential. The energy equation is derived based on Lord and Shulman theory with a single relaxation time. In the next step, by introducing dimensionless variables, the governing equations are provided in dimensionless presentation. Then these equations have been discretized using generalized differential quadrature method. Also, in order to follow the solution of the equations in time domain, Newmark method has been used. Since the system of equations is nonlinear, Picard algorithm is applied as a predictor-corrector mechanism to solve the nonlinear system of equations. Then numerical results are presented to investigate the propagation of mechanical, thermal and electric waves inside the heterogeneous sphere and also their reflection from the outer surface of the sphere. By examining the results, it can be seen that mechanical and thermal waves propagate with a limited speed, while the speed of electric wave propagation is infinite.