function nyq_rep(r,p,q)

h=p/100;
k=q/50;
x(1) = 0;
y(1) = q;
u1(1) = x(1)^3+1.6*x(1)^2+0.65*x(1)+0.05-3*x(1)*y(1)^2-1.6*y(1)^2+0.182132*exp(-r*x(1))*cos(r*y(1));
w1(1) = 3*x(1)^2*y(1)+3.2*x(1)*y(1)-y(1)^3+0.65*y(1)-0.182132*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)^3+1.6*x(i)^2+0.65*x(i)+0.05-3*x(i)*y(i)^2-1.6*y(i)^2+0.182132*exp(-r*x(i))*cos(r*y(i));
    w1(i) = 3*x(i)^2*y(i)+3.2*x(i)*y(i)-y(i)^3+0.65*y(i)-0.182132*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)^3+1.6*x(i)^2+0.65*x(i)+0.05-3*x(i)*y(i)^2-1.6*y(i)^2+0.182132*exp(-r*x(i))*cos(r*y(i));
    w2(i-101) = 3*x(i)^2*y(i)+3.2*x(i)*y(i)-y(i)^3+0.65*y(i)-0.182132*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)^3+1.6*x(i)^2+0.65*x(i)+0.05-3*x(i)*y(i)^2-1.6*y(i)^2+0.182132*exp(-r*x(i))*cos(r*y(i));
    w3(i-201) = 3*x(i)^2*y(i)+3.2*x(i)*y(i)-y(i)^3+0.65*y(i)-0.182132*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)^3+1.6*x(i)^2+0.65*x(i)+0.05-3*x(i)*y(i)^2-1.6*y(i)^2+0.182132*exp(-r*x(i))*cos(r*y(i));
    w4(i-301) = 3*x(i)^2*y(i)+3.2*x(i)*y(i)-y(i)^3+0.65*y(i)-0.182132*exp(-r*x(i))*sin(r*y(i));
end
plot(u1,w1,'b-',u2,w2,'r-',u3,w3,'g-',u4,w4,'y-');grid;