function yp=en_sust(t,y,flag,k1,k2,k3)
yp=zeros(4,1);
s=y(1);e=y(2);c=y(3);p=y(4);
yp(1)=-k1*e*s+k2*c;
yp(2)=-k1*e*s+(k2+k3)*c;
yp(3)=k1*e*s-(k2+k3)*c;
yp(4)=k3*c;