% Jacobi elliptic functions


x=-10:0.01:10;
kvec=[0,0.5,0.9,0.99,0.999,1];

%[T,E] = ellipke(k^2);

figure(1)
clf

for ip=1:6
  subplot(3,2,ip)
  k=kvec(ip);
  [Sn,Cn,Dn]=ellipj(x,k^2);
  plot(x,Cn,'LineWidth',2)
  hold on
  plot(x,cos(x),'r--',x,sech(x),'m--')
  if(ip>4)
  xlabel('x','FontSize',16); 
  end
  ylabel('cn(x;k)','FontSize',16);
  title(['k=' num2str(k)],'FontSize',16);
  set(gca,'FontSize',[16]);
end
legend('cn','cos','sech')


figure(2)
clf

for ip=1:6
  subplot(3,2,ip)
  k=kvec(ip);
  [Sn,Cn,Dn]=ellipj(x,k^2);
  plot(x,Cn.^2,'LineWidth',2)
  hold on
  plot(x,cos(x).^2,'r--',x,sech(x).^2,'m--')
  if(ip>4)
  xlabel('x','FontSize',16); 
  end
  ylabel('cn^2(x;k)','FontSize',16);
  title(['k=' num2str(k)],'FontSize',16);
  set(gca,'FontSize',[16]);
end
legend('cn^2','cos^2','sech^2')