function heumann1 = heumann1(Px) global Q f s p H1 H2 H3 h1 h2 h3 Ma P=Px(1); x=Px(2); heumann1 =[P-2*Q*f*((x/((x^2+s^2)^0.5))-((p-x)/(((p-x)^2+s^2)^0.5)))-H1-H2-H3; P*x-2*Q*f*(((x^2+s^2)^0.5)+(((p-x)^2+s^2)^0.5))-H1*(h1-x)-H2*(h2-x)-H3*(h3-x)-Ma]; end