    function z = fd(n,x,y)

    z=zeros(n,1);
    z(1)=3*(y(2)+y(1)-1/3*y(1)^3-1.3);
    z(2)=-(y(1)-0.7+0.8*y(2))/3;