function f=func3(x,p)

thetu1=x(1);
thetm1=x(2);
xu1=x(3);
delu1=x(4);

mtau1=p(1);
Dw1=p(2);
Mi=p(3);
Dku1=p(4);
Dkm1=p(5);
Dr1x=p(6);
grate=p(7);
Dxm1=p(8);
delm1=p(9);


f =[Dku1*(Mi+delu1)-thetu1/(thetu1+thetm1)*(Dr1x+xu1);
    Dkm1*(Mi/(1-mtau1)+delm1)-thetm1/(thetu1+thetm1)*(Dr1x+xu1);
    thetu1+thetm1-1+Dw1/(Dr1x+xu1+Dw1);
    xu1-(grate+delu1)*Dku1];
    

