Software to do a jacobian matrix of four variable I have to compute the jacobian wrt $\phi,\lambda_1,\lambda_2,\sigma$ of 
$C^S(\mathbf{h,x_0})=\sigma\sqrt{\mathbf{h^t\Sigma_{x_0}^{-1}h}}$
where: $$\mathbf{\Sigma_{x_0}=\Psi_{x_0}\Lambda_{x_0}\Psi_{x_0}^T } \\ \Psi_{x_0}=\begin{bmatrix}
    cos \psi(x_0)       & sin\psi(x_0)\\
    -sin\psi(x_0)       & cos \psi(x_0)\\
    \end{bmatrix}   \\ \Lambda=\begin{bmatrix}
    \lambda_1(x_0)      & 0\\
    0       & \lambda_2(x_0)\\
    \end{bmatrix} 
\ h=\begin{bmatrix}h_1 \\ h_2 \end{bmatrix}$$
Where $h$ is a known vector.
I want to compute the Jacobian wrt the four parameters, I would do through a software, I tried with wolfram alpha but it doesn't work:wolfram alpha 
 A: For typing convenience, I'll use a colon to denote the trace product, i.e.
$$A:B={\rm tr}(A^TB)$$
For the same reason, I'll use variables with Latin, rather than Greek, names 
$$\eqalign{
 L=\Lambda &= {\rm Diag}(\lambda) \cr
 P=\Psi &= I\cos\phi + K\sin\phi &\implies dP=(K\cos\phi-I\sin\phi)\,d\phi \cr
 S=\Sigma &= PLP^T \cr
 M &= S^{-1} &\implies dM=-M\,dS\,M \cr
 C^2 &= \sigma^2M:hh^T \cr
}$$
Taking the differential of the last expression
$$\eqalign{
 \frac{2C}{\sigma^2}\,dC
 &= 2M:dh\,h^T - hh^T:M\,dS\,M \cr
 &= 2Mh:dh - Mhh^TM:dS \cr
}$$
Pause at this point to note that we have found the first of the requested gradients 
$$\eqalign{
\frac{\partial C}{\partial h} &= \frac{\sigma^2}{C}\,Mh \cr
}$$
Continuing with the other half of the differential
$$\eqalign{
 \frac{2C}{\sigma^2}\,dC
 &= -Mhh^TM:dS \cr
 &= -Mhh^TM:2\,{\rm sym}(dP\,LP^T) - Mhh^TM:P\,dL\,P^T \cr
 &= -2Mhh^TMPL:dP - P^TMhh^TMP:dL \cr
 &= -2Mhh^TMPL:(K\cos\phi-I\sin\phi)\,d\phi - P^TMhh^TMP:{\rm Diag}(d\lambda) \cr
 &= -2Mhh^TMPL:(K\cos\phi-I\sin\phi)\,d\phi - {\rm diag}(P^TMhh^TMP):d\lambda \cr
}$$
Now we can write two more gradients
$$\eqalign{
\frac{dC}{d\phi}
  &= \frac{\sigma^2}{C}\,Mhh^TMPL:(I\sin\phi-K\cos\phi) \cr
\frac{\partial C}{\partial\lambda}
  &= -\frac{\sigma^2}{2C}\,{\rm diag}(P^TMhh^TMP) 
  &\implies
\frac{dC}{d\lambda_k}
 = -\frac{\sigma^2}{2C}\,e_k^TP^TMhh^TMPe_k \cr
}$$
And finally, the simplest of the gradients is
$$\eqalign{
\frac{dC}{d\sigma} &= \frac{\sigma}{C}\,h^TMh \cr
}$$
