function nyq(a,b,c,r,p,q)

h=p/100;
k=q/50;
x(1) = 0;
y(1) = q;
u1(1) = (x(1)+a)*(x(1)+b)-y(1)^2+c*exp(-r*x(1))*cos(r*y(1));
w1(1) = y(1)*(2*x(1)+a+b)-c*exp(-r*x(1))*sin(r*y(1));
for i=2:101
    x(i) = 0;
    y(i) = y(i-1)-k;
    u1(i) = (x(i)+a)*(x(i)+b)-y(i)^2+c*exp(-r*x(i))*cos(r*y(i));
    w1(i) = y(i)*(2*x(i)+a+b)-c*exp(-r*x(i))*sin(r*y(i));
end
for i=102:201
    x(i) = x(i-1)+h;
    y(i) = -q;
    u2(i-101) = (x(i)+a)*(x(i)+b)-y(i)^2+c*exp(-r*x(i))*cos(r*y(i));
    w2(i-101) = y(i)*(2*x(i)+a+b)-c*exp(-r*x(i))*sin(r*y(i));
end
for i=202:301
    x(i) = p;
    y(i) = y(i-1)+k;
    u3(i-201) = (x(i)+a)*(x(i)+b)-y(i)^2+c*exp(-r*x(i))*cos(r*y(i));
    w3(i-201) = y(i)*(2*x(i)+a+b)-c*exp(-r*x(i))*sin(r*y(i));
end
for i=302:401
    x(i) = x(i-1)-h;
    y(i) = q;
    u4(i-301) = (x(i)+a)*(x(i)+b)-y(i)^2+c*exp(-r*x(i))*cos(r*y(i));
    w4(i-301) = y(i)*(2*x(i)+a+b)-c*exp(-r*x(i))*sin(r*y(i));
end
plot(u1,w1,'b-',u2,w2,'r-',u3,w3,'g-',u4,w4,'y-');grid