What kind of software can be used to solves systems of equations? For example, I have to solve the following equations:
$$\left\{\begin{align*}
&x^2 + y^2 + z^2 = 1\\
&Ax + By + Cz = 0
\end{align*}\right.$$
for $y$ and $z$, where $x$, $A$, $B$ and $C$ are known values. I need this for production. I guess I could solve this by hand, but there are many other of those and, if did all manually, my production would stop. Is there a software that solves this kind of situation?
 A: You may consult Maple to do this for you by the following codes. However, knowing the unknown coefficients make the way of finding the variables easier.
[> solve({x^2+y^2+z^2 = 1, A*x+B*y+C*z = 0}, {y, z});

 {y = -(A*x+C*RootOf((C^2+B^2)*_Z^2+2*A*x*C*_Z+x^2*B^2+A^2*x^2-B^2))/B, z = RootOf((C^2+B^2)*_Z^2+2*A*x*C*_Z+x^2*B^2+A^2*x^2-B^2)}

A: Hot off the Mathematica presses:
> Solve[{x^2 + y^2 + z^2 - 1 == 0, A x + B y + C z == 0}, {y, z}] // Simplify

$$\begin{align}
\Bigg\{\Bigg\{y & \to -\frac{A B x+\sqrt{-C^2 \left(A^2 x^2+B^2 \left(-1+x^2\right)+C^2 \left(-1+x^2\right)\right)}}{B^2+C^2}, \\
z & \to \frac{-A C^2 x+B \sqrt{-C^2 \left(A^2 x^2+B^2 \left(-1+x^2\right)+C^2 \left(-1+x^2\right)\right)}}{C \left(B^2+C^2\right)}\Bigg\}, \\
\Bigg\{y & \to \frac{-A B x+\sqrt{-C^2 \left(A^2 x^2+B^2 \left(-1+x^2\right)+C^2 \left(-1+x^2\right)\right)}}{B^2+C^2}, \\
z & \to -\frac{A C^2 x+B \sqrt{-C^2 \left(A^2 x^2+B^2 \left(-1+x^2\right)+C^2 \left(-1+x^2\right)\right)}}{C \left(B^2+C^2\right)}\Bigg\}\Bigg\} \\
\end{align}$$
I have intentionally not simplified these equations, so as to give you a feel what Mathematica solutions are like.
A: Based on Mario Carneiro's comment, you can do just this:
Let $V = (A,B,C)$ be a vector. Your equation is asking for the set of points $(x,y,z)$ on the unit sphere that is orthogonal to $V$. 


*

*If $A = B = 0$ we have that the circle is just the circle in the $x, y$ plane. So you can parametrise as normal: $\theta \mapsto (\cos\theta, \sin\theta, 0)$ would give you points on this circle by angle. 

*If at least one of $A,B$ are not zero you can rotate the unit circle from the $x,y$ plane via an orthogonal transformation. We do this by constructing an orthogonal matrix that sends the $z$ axis to the subspace spanned by the vector $V$. Since orthogonal transformations preserve length and angles, the image of the unit circle from the $x,y$ plane under this orthogonal transformation would be the circle you want. There are lots of ways to do this. We choose one here that is algebraically simpler. 
Set the vectors $W_3 = V$, $W_2 = (B, -A, 0)$, $W_1 = (-AC,-BC,B^2+A^2)$. A direct computation will verify to you that the three are orthogonal. Set $U_3 = W_3 / \|W_3\|$ where $\|\cdot\|$ is the vector 2-norm, so $U_3 = W_3 / \sqrt{A^2 + B^2 + C^2}$. Set $U_2 = W_2 / \|W_2\|$ and set $U_1 = W_1 / \|W_1\|$. Now the matrix $O$ whose columns are given by the vectors $U_1, U_2, U_3$ is orthogonal. And it maps the $z$ axis to the subspace spanned by $U_3$, which is also the subspace spanned by $V$. 
So the following mapping would give you a parametrisation of the required circle:
$$ \theta \mapsto \cos\theta U_1 + \sin\theta U_2 $$ 
for completeness we expand the definition of $U_1$ and $U_2$ to get
$$ \theta \mapsto \cos\theta \frac{(-AC, -BC, A^2 + B^2)}{\sqrt{ (A^2 + B^2+ C^2)(A^2 + B^2)}} + \sin\theta \frac{(B,-A,0)}{\sqrt{A^2 + B^2}} $$
Note that when at least one of $A,B$ is non-zero, the denominators are non-zero and so the mapping is well defined. 
One caveat: because of the method of construction, the family of maps is discontinuous in $(A,B,C)$ as $A,B$ both approach 0. This is why we have to fill in the edge case specifically in step 1. 

In view of Dokkat's comment below, let me quickly describe the construction of $O$ and the choice of $W_1,W_2,W_3$. The first thing to remember from linear algebra is that 

If $O$ is the matrix corresponding to an orthogonal transformation, then its column vectors form an orthonormal basis of the vector space, and vice versa. 

Furthermore, if the column vectors of $O$ are $U_1, U_2, U_3$ in order, then the vector $(0,0,1)$ would be mapped to $U_3$ by the transformation $O$ by matrix multiplication. 
So to get our desired $O$ (which sends the $z$ axis to the subspace spanned by $V$), we need that $U_3$ is the unit vector corresponding to $V$. Then we just need to find $U_1, U_2$ such that $U_1, U_2, U_3$ together is an orthonormal basis. 
Given that the direction of $U_3$ is fixed, we first fix the directions of $U_1$ and $U_2$. After we fix the directions, we can just take the unit vectors in those directions to fix the length. 
To fix the directions we observe the following two facts:


*

*Suppose $A \neq 0$, the vector $(A,B,C)$ is perpendicular to the vector $(B,-A,0)$. This is an analogue of "rotating by 90 degrees in the plane" when we think about vectors in a plane. Note that we keep only the $x,y$ components and drop the $z$ component. You can do this in also higher dimensional spaces. 
We set $W_2$ to be this vector. 

*In 3D space, we have the cross product. Which has the nice property that given vectors $\vec{v}, \vec{w}$, their cross product $\vec{v}\times \vec{w}$ is perpendicular to either $\vec{v}$ or $\vec{w}$. Using this property we set $W_1 = W_2 \times W_3$. (Or maybe the negative of that, I forgot which sign I took.) 
Then the $W_1, W_2, W_3$ will be mutually orthogonal nonvanishing vectors. 
