Why coordinate axes are perpendicular to each other? Coordinate axes are chosen to be perpendicular to each other because it is convenient, calculations are easier. This is what I have read.
But in case that coordinate axes are per example 60 degrees apart could again a physical system exhibit linearity?
The superposition theorem could work?
 A: So here's the basic reason why: we like to express vectors in terms of their components, for which we invent the unit vectors $\hat e^{1,2,3}$ and components $v_{1,2,3}$ for a vector $\vec v$ such that 
$$
\vec v= v_1\hat e^1 + v_2\hat e^2 + v_3 \hat e^3.
$$
There is a bilinear combination between two vectors to form a scalar called the dot product, 
$$\vec u\cdot\vec v =\lvert \vec u\rvert\lvert \vec v \rvert\cos\theta.$$
For orthogonal unit vectors we have the nice property that 
$$\hat e^m\cdot \hat e^n=\{1\text{ if } m=n\text{ else } 0\},$$
which simplifies this expression to just
$$\vec u\cdot\vec v=u_1v_1+u_2v_2+u_3v_3.$$
This simple formula is why we use orthogonal unit vectors in practice.
There is a slightly more complicated formula in skewed  coördinate systems, or ones without unit vectors for their basis. To do this you simply have to find a new set of vectors that makes this relationship true, called the dual basis. So if your basis vectors are $\vec b^{1,2,3}$ then your dual basis is $\vec b_{1,2,3}$ such that 
 $$\vec b_{m}\cdot \vec b^n=\{1\text{ if } m=n\text{ else } 0\},$$ in other words for $\vec b^1$, say, you look at the plane spanned by $\vec b^{2,3}$, identify a vector normal to it $\vec n$, and then rescale it until its dot product with $\vec b^1$ is $1$, and that then is $\vec b_1$. So the dual basis for some periodic lattice of points tells you about the planes of that lattice, via their normal vectors.
Then every vector gets two sets of components, the regular and the dual components,
$$\begin{align}
\vec v&= v_1\vec b^1 + v_2\vec b^2 + v_3 \vec b^3
&=v^1\vec b_1 + v^2\vec b_2 + v^3 \vec b_3,
\end{align}$$
at which point the earlier easy property is restored as long as we always pair a lower index with an upper index,$$
\vec u=u_1v^1+u_2v^2+u_3v^3=u^1v_1+u^2v_2+u^3v_3.$$
A: 
The   vector $\vec{R}$ can described with the  vectors basis $\vec{q}_1\,,\vec{q}_2$ :
$$\vec{R}_s= a_1\begin{bmatrix}
     1\\
     0\\
  \end{bmatrix}+a_2\begin{bmatrix}
     \cos(\theta)\\
     \sin(\theta)\\
  \end{bmatrix}=a_1\,\vec{q}_1+a_2\,\vec{q}_2
\tag 1$$
or with the  vectors basis $\vec{e}_1\,,\vec{e}_2$
$$\vec{R}_c= c_1\begin{bmatrix}
     1\\
     0\\
  \end{bmatrix}+c_2\begin{bmatrix}
     0\\
     1\\
  \end{bmatrix}=c_1\,\vec{e}_1+c_2\,\vec{e}_2
\tag 2$$
both case the superposition theorem is valid
Edit 
what's about the EOMs?
I)  generalized coordinates are $x$ and $y$
with $c_1=x$ and $c_2=y$ we obtain for the position Vector $\vec{R}_c$ equation (2)
$$\vec{R}_c= x\begin{bmatrix}
     1\\
     0\\
  \end{bmatrix}+y\begin{bmatrix}
     0\\
     1\\
  \end{bmatrix}
$$ 
Kinetic energy 
$T=m\,\frac{1}{2}\dot{R}_c^T\,\dot{R}_c$ and
Potential energy
$V=-m\,g\,y$
$\Rightarrow\quad$  EOMs
$$ \begin{bmatrix}
     \ddot{x}\\
     \ddot{y}\\
  \end{bmatrix}=\left[ \begin {array}{c} 0\\ -g\end {array}
 \right] 
\tag 3$$
II)  generalized coordinates are $q_1$ and $q_2$
With $a_1=q_1$ and $a_2=q_2$ we obtain for the Position Vector $\vec{R}_s$ equation (1)
$$\vec{R}_s= q_1\begin{bmatrix}
     1\\
     0\\
  \end{bmatrix}+q_2\begin{bmatrix}
     \cos(\theta)\\
     \sin(\theta)\\
  \end{bmatrix}$$
Kinetic energy 
$T=m\,\frac{1}{2}\dot{R}_s^T\,\dot{R}_s$ and
Potential energy
$V=-m\,g\,q_2\,\sin(\theta)$
$\Rightarrow\quad$ EOMs
$$ \begin{bmatrix}
     \ddot{q}_1\\
     \ddot{q}_2\\
  \end{bmatrix}=\left[ \begin {array}{c} {\frac {g\cos \left( \theta \right) }{\sin
 \left( \theta \right) }}\\ -{\frac {g}{\sin \left( 
\theta \right) }}\end {array} \right] \tag 4
$$
The generalized accelerations equation (3) and (4) are not equal. we can  transfer equation (4) to get the accelerations in the orthogonal coordinates  $\vec{e}_1\,,\vec{e}_2$
$$\begin{bmatrix}
     \ddot{x}\\
     \ddot{y}\\
  \end{bmatrix}=\left[ \begin {array}{cc} 1&\cos \left( \theta \right) 
\\ 0&\sin \left( \theta \right) \end {array}
 \right]\,\begin{bmatrix}
     \ddot{q}_1\\
     \ddot{q}_2\\
  \end{bmatrix}=\left[ \begin {array}{c} 0\\ -g\end {array}
 \right] 
\tag 5$$
so the accelerations equation (3) and (5) are now the same.
Conclusion
The equations of motion described with  skew vectors basis compare to the  equations of motion  described with orthogonal vectors basis are not the same!.
