I have a spline generation system using a cubic Bézier algorithm.
I created a tool that splits the spline and adds a new point where the user clicked on the spline, but it causes the spline to deform from its current shape.
How can I calculate Bézier control point positions in 3D space to avoid deformation, as software like 3ds Max and Blender do?
What you're looking for is called "de Casteljau's algorithm"
First, you need to find the segment and the parameter value \$t\$ at the point that was clicked. Search "closest point on cubic Bézier spline" for algorithms to find that.
Now you have a segment with control points \$\{ P_1, P_2, P_3, P_4 \}\$ and a parameter value \$0 < t < 1\$ (if \$t\$ is exactly zero or one, then control point \$P_1\$ or \$P_4\$ respectively already split the spline at the clicked point), at which the Bézier segment passes through the clicked point \$P_t\$.
We can now form new points by interpolating between the existing control points by a factor of \$t\$:
$$ P_5 = (1 - t) P_1 + t P_2\\ P_6 = (1 - t) P_2 + t P_3\\ P_7 = (1 - t) P_3 + t P_4\\ $$
and again:
$$ P_8 = (1 - t) P_5 + t P_6\\ P_9 = (1 - t) P_6 + t P_7\\ $$
(If we wanted, we could also find \$P_t = (1 - t) P_8 + t P_9\$ this way, but by this stage you probably already know it)
Your two new Bézier segments are:
$$ \{P_1, P_5, P_8, P_t\}\\ \{P_t, P_9, P_7, P_4\} $$
