## 1 Introduction  
Demand for highly effective, efficient, and reliable numerical methods has grown increasingly high for solving option-trading modeling equations involving cross-derivative terms. However, desirable computational procedures are, in general, difficult to obtain, due to challenges from the participation of cross-derivatives [1; 2]. This motivates our study. In this investigation, targeted at European options that can only be exercised on dates of maturity, we propose and analyze a new and dynamically balanced up-downwind finite difference method in the pursuit. In the early 1970s, [PERSON], [PERSON], and [PERSON] introduced the popular Black-Scholes-Merton (BSM) model [3; 4]. Under their consideration, stock prices were assumed to follow geometric Brownian motion, while the volatility of the stock prices was fixed and no sudden jumps occurred. However, classic BSM models often cannot fit ideally into the market data observed nowadays [4]. This may be due to the fact that, in modern financial markets, not only are stock prices subject to risk, but also the estimate of riskiness is typically subject to significant uncertainty. The Black-Scholes model does not adequately take into account essential characteristics of market dynamics, such as fat tails, skewness of the distribution of log returns, and the correlation between the value of the underlying and its volatility. It has also been observed that the volatility starts to fluctuate when the market reacts to new information [5]. To incorporate an additional source of randomness into an option pricing model, [PERSON] proposed a more refined approach, based on the concept of stochastic volatility [5; 6; 7]. A closed-form solution of the model was also obtained by [PERSON], under a set of specific boundary and initial values for assets of the European type [6]. However, to meet a growing demand for American options and other assets, pricing equations often need to be placed together with more realistic initial boundary conditions, or even free boundary conditions. Closed forms of solutions are, in general, unavailable. Thus, numerical approximations of such solutions have become important and necessary. This paper is concerned with European options. The scheme we developed here can be extended directly to price American options. However, due to the free boundary conditions associated with the American options, the stability and convergence analysis of the scheme becomes rather complicated. We are still working on the theoretical results behind applying our scheme to American options. There have been numerous recent publications on the numerical solution of Heston modeling equations. For instance, certain first-order up-downwind algorithms were proposed and studied by [PERSON] and [PERSON] [2]. Stability analysis was carried out via standard von Neumann analysis for Cauchy problems or problems with periodic boundary conditions [8; 9]. Although numerical stabilities have been under investigation for even high order schemes on nonuniform grids [4], rigorous analysis is only available in cases where the cross-derivative terms are neglected. The challenge for stability analysis persists, whenever a cross-derivative or more general boundary data structure exists [2; 10]. However, cross-derivatives are essential to partial differential equations modeling a Heston Process. Further, Heston modeling formulations also require more realistic Dirichlet, Neumann, or mixed boundary conditions [6; 11]. These have motivated our approaches. In this paper, we are particularly interested in computations based on a Heston put option model [4; 8; 9; 12; 13; 14]. We are primarily interested in a linearly stable finite difference method, based on nonuniform grids. Our intention is to effectively reduce the computational costs and raise the algorithmic efficiency, by way of application of appropriate adaptive mechanisms. These results can be extended for call options, in similar ways.