function bar4()
clear; close all

a1=10;a2=5;a3=7;a4=2;
psi0=-30; % change for different zeros
A=[0,0];
B=[a1,0];
hold on
axis equal
for i = 1:60
phi(i)=i*pi/30;
psi(i) = fzero(@bar4,psi0,[],phi(i),a1,a2,a3,a4);
C=[a1+a2*cos(psi(i)),a2*sin(psi(i))];
D=[a4*cos(phi(i)),a4*sin(phi(i))];
XX=[A;B;C;D;A];
X1 = XX(:,1);Y1 = XX(:,2);
plot(X1,Y1,'k-')
x = (a1+a2*cos(psi)+a4*cos(phi))/2;
y = (a2*sin(psi)+a4*sin(phi))/2;
plot(x,y,'ro')
pause(.1)
plot(XX(:,1),XX(:,2),'y-')
end
plot(XX(:,1),XX(:,2),'k-')
hold off

function z = bar4(psi,phi,a1,a2,a3,a4)
z=a1^2+a2^2-a3^2+a4^2-2*a2*a4*cos(phi-psi)+...
2*a1*(a2*cos(psi)-a4*cos(phi));