weird branching with elliptic12i (Matlab/Octave)

[dezanche]

I'm using your library in Octave to compute conformal maps around rectangular domains, but I get weird results like this instead of smooth streamlines:

L=10;      % number of horizontal streamlines
p=0.4;    % parameter that determines rectangle's aspect ratio
[Kp Ep] = ellipke(1-p); % K' and E'
[K E] = ellipke(p); % K and E
% dimensions of rectangle (eq. 1.13.13 Ivanov, Trubetskov)
rect_Y = (Ep-p*Kp)./sqrt(p); 
rect_X = (E-(1-p)*K)./sqrt(p);
% Im>0 half plane
[Re, Im] = meshgrid(0.2*[-7:7],0.01+0.1*[0:L]);
w = Re + 1i*Im;
% calculation of E(t,1/p):
[F, Ei, Z] = elliptic12i(asin(w)/sqrt(p), p); 
E1p = Ei/sqrt(p)-(1/p-1)*asin(w);    % Abramowitz & Stegun eq. 17.4.16
figure(1);
subplot(1,2,1)
plot(w(1,:),'k.');
hold on
plot(w(1,:),'k-');
for l = 2:L
  plot(w(l,:),'k.');
  plot(w(l,:),'k-');
end
hold off
axis equal;
title('Cartesian grid');
subplot(1,2,2)
plot(E1p(1,:),'k.');
hold on
plot(E1p(1,:),'k-');
plot(rect_X*[-1 -1 1 1],rect_Y*[-1 0 0 -1],'k--');
plot(rect_X*[-3 3],rect_Y*[-1 -1],'k--');
for l = 2:L
  plot(E1p(l,:),'k.');
  plot(E1p(l,:),'k-');
end
hold off
axis equal;
title('transformed grid');

Any ideas on how to fix this?
I checked the other functions in my code and they are not introducing these discontinuities. I suspect it has something to do with how elliptic12i assembles the output from the real results of elliptic12.