I have coordinates of triangle's points (x1, y1, x2, y2, x3, y3) and line equation y = kx + b.
My idea is to convert line equation to Ax + By + C and find line equation for each side of triangle (A1x + B1y + C1, A2x + B2y + C2, A3x + B3y + C3). Then simply solve the system of equations for each side of triangle and the line. But I'm not sure how to check that the points of intersection are on the triangle sides.
What is the right way to find points where line intersect the triangle?
Besides finding these intersections, you also have to check that they lie between the corresponding vertices. I don’t know if this is the “right” way, but a relatively simple way to perform this check in the course of finding the intersections is to do the calculations in barycentric coordinates.
On the extended Euclidean plane (the projective plane with some extra geometry) we can represent the line $y=mx+b$ by the homogeneous coordinate row vector $\mathbf l = (m,-y,b)$. Since the barycentric coordinate system of the triangle is another homogeneous coordinate system, we can treat it as another projective plane, so that the Cartesian-to-barycentric coordinate transformation is just a homography between the two planes. Letting $$M=\begin{bmatrix}x_1&x_2&x_3\\y_1&y_2&y_3\\1&1&1\end{bmatrix}$$ the line’s representation in barycentric coordinates is its image $\mathbf l'=\mathbf lM$ under this homography. Now, the barycentric equations of the extensions of the triangle’s sides are simply $\lambda_i=0$, so the barycentric coordinates of the intersection with $\mathbf l'$ are its cross products with the standard unit vectors. That is, if $\mathbf l'=(\mu,\nu,\tau)$, then the (unnormalized) barycentric coordinates of its intersections with the side extensions are $$\begin{bmatrix}0\\-\tau\\\nu\end{bmatrix},\begin{bmatrix}\tau\\0\\-\mu\end{bmatrix},\begin{bmatrix}-\nu\\\mu\\0\end{bmatrix}.$$ (I’m using column vectors for points to distinguish them from the covariant vectors that represent lines.) Compare the signs of the nonzero components of each point: if they are the same, then that point lies on the triangle. If you were to normalize the coordinates, they would both be positive. (Note that if the coordinates sum to $0$, then the line is parallel to that side, but you’ll have already filtered that case out because the signs are different.) If two of the point’s barycentric coordinates are zero, then the line passes through a vertex. If all three are zero, the line coincides with the corresponding side of the triangle.
To convert the barycentric coordinates $\mathbf\lambda$ of a point back to Cartesian, compute $M\mathbf\lambda$ and dehomogenize by dividing through by the third component.†
For example, with the vertices $(0,0)$, $(1,0)$ and $(2,3)$ and the line $y=1-x/3$, we have $$\mathbf l M = \begin{bmatrix}-\frac13&-1,1\end{bmatrix} \begin{bmatrix}0&1&2\\0&0&3\\1&1&1\end{bmatrix} = \begin{bmatrix}1&\frac23&-\frac83\end{bmatrix},$$ so that the three intersection points are $\left(0,-\frac83,-\frac23\right)^T$, $\left(\frac83,0,1\right)^T$ and $\left(\frac23,-1,0\right)^T$. The signs of the nonzero components of the third point don’t match, so that point doesn’t lie on the triangle; the other two do. Converting the first one back to Cartesian coordinates, we have $$M\begin{bmatrix}0\\-\frac83\\-\frac23\end{bmatrix} = \begin{bmatrix}4\\2\\\frac{10}3\end{bmatrix},$$ which dehomogenizes to $\left(\frac65,\frac35\right)$. Similarly, the third point maps to the Cartesian coordinates $(3,0)$, which is obviously not on this triangle.
† Equivalently, you can of course normalize the barycentric coordinates first and then compute $\lambda_1\mathbf p_1+\lambda_2\mathbf p_2+\lambda_3\mathbf p_3$, where the $\mathbf p_i$ are the Cartesian coordinates of the triangle’s vertices. Effectively, you’d be performing linear interpolations between pairs of vertices.