  // exact solution of the
  // nonlinear BVP
  // y''=3/2 y^2, y(0)=4, y(1)=1
  
  j=2;
  h=10^(-j);
  x=0:h:1;

  // first solution
  y1=4.0./(1+x).^2;

  // second solution
  k=sqrt(2.0+sqrt(3.0))/2.0;
  m=k^2;
  C1=4.30310990;
  C2=2.33464196;
  // cn=real(%cn(C1*x-C2,m));
  // function cn does not exist in Scilab, therefore:
  cn=sqrt(1-real(%sn(C1*x-C2,m)).^2); 
  y2=C1^2*((1.0-cn)./(1.0+cn)-1.0/sqrt(3.0));

  xbasc()
  plot([y1;y2]) 
