What I need is to generate an SVG file while having a series of (x,y) ready.
P0(x0,y0)
P1(x1,y1)
P2(x2,y2)
P3(x3,y3)
P4(x4,y4)
P5(x5,y5)
...
I need to make a Bézier curve from them so I need to calculate control points (mid-points) for them. This image explains what I am exactly looking for:

I have points P0, P1, P2 and P3 ready. What I need is to calculate control points C1 and C2. The curve does not pass them. But it bends toward them.
I need a formula which gives me C1 and C2 in clear direct form:
C1= fomula1 (P0,P1,P2,P3)
C2= fomula2 (P0,P1,P2,P3)
I was thinking about least square method and some other methods but I had no idea how to implement them.
Your problem, as stated, does not have a unique solution. Suppose that point $P_j$ is at location $(3j, 0)$, for each integer $j$, so that they're equi-spaced on the $x$-axis. Now let $y$ be any real number. Then by adding control points at locations $$(6i+1, y)\\ (6i+2, y)\\ (6i + 4, -y)\\ (6i+5, -y)$$ for each integer $i$, you get two "control points" between any two of your original points. For instance, near the origin, for $y = 2$, we have points $$ (-3, 0) \leftarrow (\mathrm{one~ of~ the}~ P_i)\\ (-2, -2)\\ (-1, -2)\\ (0, 0) \leftarrow (\mathrm{one~ of~ the}~ P_i)\\ (1, 2)\\ (2, 2)\\ (3, 0) \leftarrow (\mathrm{one~ of~ the}~ P_i) $$ These determine two Bezier segments that glue up nicely at the origin, with a slope of $2$ at the origin.
You may say "But it's obvious that the control points in this case should be on the $x$-axis!" and I say "but your problem statement doesn't require it." Indeed, I chose this example because it was easy to write, but given any set of $P_i$, I can again find an infinitely family of ways to place the intermediate control points so as to join the $P_i$ with Bezier segments.
I'm going to suggest that you consider looking at Catmull-Rom splines, which are piecewise cubics passing through a sequence of points like your $P_i$. Each segment of a CR-spline can be expressed as a Bezier curve, because the Bezier basis functions span the space of cubic curves. One detailed reference on this is Computer Graphics: Principles and Practice, 3rd edition, of which I am a coauthor, but there are plenty of other references as well.
Here are somewhat brief details on CR spline construction from a sequence of points $M_0, M_1, \ldots$. I'm going to describe how to find the control points for the part of the curve between $M_1$ and $M_2$, so as to avoid any negative indices. The four control points will be $P_0, P_1, P_2, P_3$. Two of these are easy: \begin{align} P_0 &= M_1 \\ P_3 &= M_2 \end{align} so that the Bezier curve starts and ends at $M_1$ and $M_2$, respectively.
The other two are only slightly trickier. We compute \begin{align} v_1 &= \frac{1}{2} (M_2 - M_0)\\ v_2 &= \frac{1}{2} (M_3 - M_1) \end{align} which are the velocity vectors at $M_1$ and $M_2$. We then have \begin{align} P_1 &= P_0 + \frac{1}{3} v_1 = M_1 + \frac{1}{6} (M_2 - M_0)\\ P_2 &= P_4 - \frac{1}{3} v_2 = M_2 - \frac{1}{6} (M_3 - M_1) \end{align}
Applying these rules in the example I gave earlier, with $$ M_i = (3i, 0) $$ we have \begin{align} M_0 &= (0, 0)\\ M_1 &= (3, 0)\\ M_2 &= (6, 0)\\ M_3 &= (9, 0) \end{align} so that \begin{align} P_0 &= (3, 0)\\ P_3 &= (6, 0)\\ v_1 &= \frac{1}{2}((6, 0) - (0, 0)) = (3, 0)\\ v_2 &= \frac{1}{2}((9, 0) - (3, 0)) = (3, 0)\\ P_1 &= P_0 + \frac{1}{3} v_1 = (3, 0) + (1, 0) = (4, 0)\\ P_2 &= P_3 - \frac{1}{3} v_2 = (6, 0) - (1, 0) = (5, 0) \end{align} as expected.
Hint for the start and end points: Assuming you have a sequence of points $$ M_0, M_1, \ldots, M_{n-1} $$ you can let \begin{align} M_{-1} &= M_0 - (M_1 - M_0) = 2M_0 - M_1 \\ M_n &= M_{n-1} + (M_{n-1} - M_{n-2}) = 2M_{n-1} - M_{n-2} \end{align} to extend your list just enough that the CR scheme above provides interpolation all the way from $M_0$ to $M_{n-1}$.