algorithm to calculate the control points of a cubic Bezier curve

I have all points where my curve pass through, but I need to get the coordinates of the control points to be able to draw the curve. How can I do to calculate this points?

-

When what you already have is a set of points where your curve must pass through, Bézier is not what you want; you should be using a parametric (cubic) spline.

Assuming that your points $\mathbf{P}_i=(x_i,y_i)$, $i=1\dots n$ represent a general curve (as opposed to a function $f(x)$, for which simply applying the usual cubic spline algorithm on your points suffices), Eugene Lee proposed centripetal parametrization to generate suitable parameter values $t_i$ associated with your $\mathbf{P}_i$. The prescription for generating the $t_i$ (in its most general form) is

where $\left\| \cdot \right\|$ is the (Euclidean) length, $e$ is an adjustable exponent in the interval $[0,1]$, and $t_1=0$. (A usual value for $e$ is 0.5, but $e=1$ is sometimes used as well.)

From this, one applies the usual cubic spline algorithm to the sets $(t_i,x_i)$ and $(t_i,y_i)$, from which you now have your parametric spline. (The periodic spline is recommended for closed curves, and the "not-a-knot" spline for all other cases.)

The MATLAB Spline Toolbox has support for this parametrization scheme, though it shouldn't be too hard to write your own implementation.

-
If you have only 4 points, then you will need only one Bezier cubic segment. Suppose your known points (through which you want the curve to pass) are $Q_0, Q_1, Q_2, Q_3$. First you have to choose 4 parameter values $t_0,t_1,t_2,t_3$ to assign to these points. The centripedal approach described in the other answer is good. In simple situations, where your points are fairly equidistant, a much simpler choice is just $(t_0,t_1,t_2,t_3) = (0, 1/3, 2/3, 1)$. Then you have to solve a system of 4 linear equations, as follows. Suppose $P_0,P_1, P_2, P_3$ are the (unknown) control points of the Bezier curve, and denote the curve by $C(t)$. Then we want $C(t_i) = Q_i$ for $i=0,1,2,3$. This means $$\sum_{j=0}^3{b_j(t_i)P_j} = Q_i$$ for $i=0,1,2,3$, where $b_0, b_1,b_2,b_3$ are the Bernstein polynomials of degree 3. Solve these equations to get $P_0,P_1, P_2, P_3$. If you always use $(t_0,t_1,t_2,t_3) = (0, 1/3, 2/3, 1)$, then the coefficient matrix of the linear system is fixed, and you can just invert it once (exactly), and save the answer. This gives you a simple formula relating the $P$'s and $Q$'s.
If you have more than 4 points, then you can use a Bezier curve with degree higher than 4. The calculations are analogous to those shown above. But a better approach is to use a "spline" that consists of several Bezier cubic segments joined end-to-end. There are many ways to compute splines. One of the easiest (and the most popular amongst graphics folks) is the Catmull-Rom spline. It gives you simple formulae for the control points in terms of the given points $Q_i$.