M=[1,0,0; 0,1,0; 0,0,0]; options=odeset; options.Mass=M;
x0=[0.8; 0.1; 0.1]; [t,x]=ode15s(@c3eqdae,[0,20],x0,options);
plot(t,x)

x0=[0.8; 0.1]; [t1,x1]=ode45(@c3eqdae1,[0,20],x0);
plot(t1,x1,t1,1-sum(x1'))