tricl(theta,y)

function yp = tricl(theta,y)

k = [ 240 30 1];
p0 = 2;
v = 1;
p = p0 -y(2)-2*y(3)...
    -3*(1-y(1)-y(2)-y(3));
q = y(1);
r = y(2);
s = y(3);
yp = [ -k(1).*p.*q;
       (k(1).*q - k(2).*r).*p;
       (k(2).*r -k(3).*s).*p;
     ]./v;