a=[-8/3 0 0; 0 -10 10; 0 28 -1];
b=[-8/3 0 0; 0 -10 10; 0 28 -1];
c=[-8/3 0 0; 0 -10 10; 0 28 -1];
d=[-8/3 0 0; 0 -10 10; 0 28 -1];
y=[20 5 6]';
z=[(20+.000001),(5+.000001),(6+.000001)]';
z1=[(20+.000002),(5+.000002),(6+.000002)]';
z2=[(20+.000003),(5+.000003),(6+.000003)]';
h=.01;
p=plot3(y(1),y(2),y(3),'.c', ...
        'EraseMode','background','MarkerSize',10)
axis([0 50 -25 25 -25 25])
hold on
q=plot3(z(1),z(2),z(3),'.r', ...
        'EraseMode','background','MarkerSize',10)
r=plot3(z1(1),z1(2),z1(3),'.g', ...
        'EraseMode','background','MarkerSize',10)
s=plot3(z2(1),z2(2),z2(3),'.b', ...
        'EraseMode','background','MarkerSize',10)
while 1
   a(1,3)=y(2);
   a(3,1)=-y(2);
   b(1,3)=z(2);
   b(3,1)=-z(2);
   c(1,3)=z1(2);
   c(3,1)=-z1(2);
   d(1,3)=z2(2);
   d(3,1)=-z2(2);
   ydot=a*y;
   zdot=b*z;
   z1dot=c*z1;
   z2dot=d*z2;
   y=y+h*ydot;
   z=z+h*zdot;
   z1=z1+h*z1dot;
   z2=z2+h*z2dot;
   set(p,'XData',y(1),'YData',y(2),'ZData',y(3))
   set(q,'XData',z(1),'YData',z(2),'ZData',z(3))
   set(r,'XData',z1(1),'YData',z1(2),'ZData',z1(3))
   set(s,'XData',z2(1),'YData',z2(2),'ZData',z2(3))
   drawnow
end