% Jacobi elliptic functions

set(groot,'defaulttextinterpreter','latex');
set(groot,'defaultAxesTickLabelInterpreter','latex');
set(groot,'defaultLegendInterpreter','latex');
set(groot,'defaultAxesFontSize',14);

x=linspace(-10,10,1000);
kvec=[0,0.5,0.9,0.99,0.999,1];
mvec=kvec.^2;

%[T,E] = ellipke(mvec);

figure(1)
clf

for ip=1:6
  subplot(3,2,ip)
  m=mvec(ip);
  k=kvec(ip);
  [Sn,Cn,Dn]=ellipj(x,m^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)
  m=mvec(ip);
  k=kvec(ip);
  [Sn,Cn,Dn]=ellipj(x,m);
  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$')