Is there a general systematic approach how to calculate the individual terms in an operator product expansion? Is there a general systematic procedure or approach to obtain the analytic functions $c_{ijk}$
as well as the corresponding operators $$A_i(z)$$ that appear in the operator product expansion (OPE)
$$
A_i(z)A_j(w) = \sum\limits_k c_{ijk}(z-w)A_k(w)
$$
for any given two operators $A_i(z)$ and $A_i(z)$? 
I have followd (a limited number of) examples of OPE calculations that involved things like expanding until some correlators for which the expression is known appear and power series expansions of derivatives of operators, but from this I was not able to discern the general way to go. 
 A: The available methods depend on what you know about the CFT (or another quantum field theory).
For example, there are nice classes of two-dimensional conformal field theories (I mean "minimal models" etc.) that are almost uniquely (up to several integer-valued labels) determined by the conformal symmetry. All the dynamics – correlators etc. – of these CFTs are fully encoded by the OPEs which means that all this dynamics is fully encoded in the coefficients $c_{ijk}$. These coefficients may be calculated as solutions to the constraints equivalent to the conformal symmetry. Once you calculate them, you should view them as a major part of the definition of the CFT – they're the most fundamental numbers defining the identity of the CFT so it's futile to try to calculate them from something more fundamental.
If you have more constructive CFTs, the CFTs may be essentially "free fermions" or "free bosons". In such a case, the OPEs may be calculated by "Wick contractions" of some kind. For more general (but non-free) CFTs, it's helpful to notice that the operators $A_j(w)$ are in one-to-one correspondence to the states in the CFT (quantized on a circular space times infinite time) and the structure coefficients $c_{ijk}$ may also be calculated from the action of the operator $A_i$ on the state corresponding to the equally named operator $|A_w\rangle$.
A: Lubo$\check{\mbox{s}}$ has already given the answer. However, I would like to share whatever I personally understand about notion of OPE.
First of all, I think, one should differentiate between notion of "operator product" and "operator product expansion". It is similar to difference between cross product of two vectors, and expansion of this cross product in a basis. Secondly this notion should not be confused to be specific only to 2d conformal field theories (though it seems to be of some practical value only in these theories).
In Hilbert space set up, operator product $O_{12}(x,y)$ of two quantum fields $O_1(x)$ and $O_2(y)$ is defined to be the time ordered product  $O_{12}(x,y)=T(O_1(x)O_2(y))$ for $x\neq y $. If we know how $O_1(x)$, and $O_2(y)$ act on given state space then in principle one can define $T(O_1(x)O_2(y))$ as an operator and so in such situation operator product is not difficult to compute. In case of some two dimensional conformal field theories singular part of the operator product of two quantum fields has the same information as contained in algebra of their modes and so usually only the singular part is mentioned in expression for the operator product of given fields. See references [1], and [2] for treatment of 2d CFT (and in particular for introduction to notion of OPE) in Hilbert space formalism.  Notion of operator product expansion can also be taken as a starting point for defining 2d CFT's (see vertex operator algebra). 
In path integral set up, operator product of two fields $O_1(x)$ and $O_2(y)$ is defined to be that field $O_{12}(x,y)$ on $M\times M-diagonal$ ($M$ being the spacetime) which is possibly singular along the diagonal and which is such that for any other fields $\phi_1 (x_1)$ ,...., $\phi_n(x_n)$ (with points $x_1$, ..., $x_n$ all distinct and away from $x$ and $y$) the correlation function of the product of fields
$((O_1(x)O_2(y) - O_{12}(x,y))  \phi_1 (x_1) .... \phi_n(x_n)$ 
goes to zero when $x\rightarrow y$.
So as follows from the definition, in path integral formalism, operator product of two fields is defined only within a regular terms which goes to zero as $x\rightarrow y$. Also in order to know operator product of two fields one needs to compute behavior of the product of given fields in all possible correlation functions. This seems to be impossible, unless given theory is simple in some sense. I personally know of only two such cases where its practically possible to compute (at least singular part of) operator product of some of the fields :- 
Case i) Given theory is free :
For a free theory operator product of "basic" dynamical fields is given by their Green's function (~ time ordered vacuum expectation value). Operator product of normal ordered fields can be computed using Wick's theorem.
Case ii) There is enough symmetry that we can make strong general statements about form of correlation functions. 
This criteria of presence of enough symmetry is e.g. met in two dimensional conformal field theories. In these theories it is possible to compute some operator products by making use of ward identities corresponding to the symmetry of the theory. Roughly speaking, Ward identities are differential equations satisfied by all correlation functions involving a particular set of fields (generators of underlying symmetry of the theory) and so can be solved to get some operator products. See reference [3] for a derivation of operator products involving energy momentum tensor and primary fields from some assumption regarding behavior of correlation functions under conformal symmetry. 

[1] Jurgen A. Fuchs, Affine Lie algebras and quantum groups, Chapter 3. 
[2] Jurgen A. Fuchs, Lectures on conformal field theory and Kac Moody algebras.
[3] K. Gawedzki, Conformal field theory, chapter (2) in "Quantum fields and strings : A course for mathematicians, Vol 2".
