%I'm going to make sure that I'm calculating the Jacobians correctly 

%The S matrices
S11=sym('1');           S12=sym('0');           S13=sym('0');
S21=sym('-u/p');        S22=sym('1/p');         S23=sym('0');
S31=sym('a*b');          S32=sym('-u*b');         S33=sym('b');

Si11=sym('1');          Si12=sym('0');          Si13=sym('0');
Si21=sym('u');          Si22=sym('p');          Si23=sym('0');
Si31=sym('a');          Si32=sym('p*u');         Si33=sym('1/b');

%The C matrices
C11=sym('1');           C12=sym('0');           C13=sym('-1/c^2');
C21=sym('0');           C22=sym('p*c');          C23=sym('1');
C31=sym('0');           C32=sym('-p*c');         C33=sym('1');

Ci11=sym('1');          Ci12=sym('1/(2*c^2)');  Ci13=sym('1/(2*c^2)');
Ci21=sym('0');          Ci22=sym('1/(2*p*c)');  Ci23=sym('-1/(2*p*c)');
Ci31=sym('0');          Ci32=sym('1/2');        Ci33=sym('1/2');

%Eigenvalue matrix      diagonal
L11=sym('u');           L12=sym('0');           L13=sym('0');
L21=sym('0');           L22=sym('u+c');         L23=sym('0');
L31=sym('0');           L32=sym('0');           L33=sym('u-c');


S=[S11 S12 S13; S21 S22 S23; S31 S32 S33];
Si=[Si11 Si12 Si13; Si21 Si22 Si23; Si31 Si32 Si33];
C=[C11 C12 C13; C21 C22 C23; C31 C32 C33];
Ci=[Ci11 Ci12 Ci13; Ci21 Ci22 Ci23; Ci31 Ci32 Ci33];
L=[L11 L12 L13; L21 L22 L23; L31 L32 L33];


Eigvecmat=simplify(C*S);
Eigvecmatin=simplify(Si*Ci);

A=expand(Eigvecmatin*L*Eigvecmat)
A
%A=Si*Ci*L*C*S;
%b=simplify(A);
%simple(b)