% here are the Pauli spin matrices and the Dirac gamma matrices


z=zeros(2);
I4=eye(4);

% Define the Pauli spin matrices

s0=eye(2);
s1=[0 1;1 0];
s2=[0 -i; i 0];
s3=[1 0; 0 -1];

% Define the gamma matrices

gd0=[s0 z  ;z -s0  ];
gd1=[z s1  ; -s1 z ];
gd2=[z s2  ; -s2 z ];
gd3=[z s3  ; -s3 z ];
gd5=i*gd0*gd1*gd2*gd3;

% The next 6 lines show anti-commutativity of the first 4 gammas

gd0*gd1+gd1*gd0
gd0*gd2+gd2*gd0
gd0*gd3+gd3*gd0
gd1*gd2+gd2*gd1
gd1*gd3+gd3*gd1
gd2*gd3+gd3*gd2

% Now define the rest of the 16 proposed basis matrices

g01=gd0*gd1;
g02=gd0*gd2;
g03=gd0*gd3;
g12=gd1*gd2;
g13=gd1*gd3;
g23=gd2*gd3;

g012=gd0*gd1*gd2;
g013=gd0*gd1*gd3;
g023=gd0*gd2*gd3;
g123=gd1*gd2*gd3;

g0123=gd0*gd1*gd2*gd3;

% We are now going to perform symbolic manipulations.
% In octave you must have the symbolic package installed and 
% then invoke it with "pkg load symbolic" 
% In MATLAB you must have the symbolic math toolbox installed

% Create the symbolic entities.

syms a b c d
syms x0 x1 x2 x3 x4 x5 x6 x7 x8 x9 x10 x11 x12 x13 x14 x15

Y=a*gd0+b*gd1+c*gd2+d*gd3

simplify(Y*Y)

X=x0*gd0+x1*gd1+x2*gd2+x3*gd3+x4*g01+x5*g02+x6*g03+x7*g12+x8*g13+x9*g23...
+x10*g012+x11*g013+x12*g023+x13*g123+x14*g0123+x15*I4

% To solve the last exercise re-run the code above with 
% redefined gd1,gd2 and gd3 as below.

gd1=i*gd1;
gd2=i*gd2;
gd3=i*gd3;