We have applied a method of characteristics to a large-scale computational problem of strongly time-dependent, two-dimensional, incompressible, infinite Prandtl number thermal convection for a non-Newtonian power law rheology. Bicubic splines are used for the spatial discretization for both temperature and stream function fields. A predictor-corrector type scheme with the Lagrangian formulation of the total time derivative in the energy equation allows one to model flows with advection strongly dominating over diffusion. The non-Newtonian rheological law introduces nonlinearity into the momentum equation. Therefore, iterations are necessary to obtain the solution at each time step. This aspect is responsible for the very CPU intensive nature of this problem. Extrapolation of the stream function fields from the previous time steps gives a good initial guess. Thus only a few Picard iterations are usually necessary to obtain solution of the momentum equation. We have observed the transition from mildly time-dependent to the strongly chaotic or turbulent non-Newtonian convection. The regular cellular structure of convection is broken in turbulent regime. The thermal plumes become disconnected and the interior of the convecting volume is filled with the detached blobs of hot and cold material. Non-Newtonian convection is far more chaotic than Newtonian carrying away the same amount of heat from a box. Rising non-Newtonian plumes exhibit much greater curvature in their ascent than Newtonian ones and are strongly attracted by descending currents at the top. Increasing the power-law index sharpens the chaotic behavior of the flow.