# This version works!

#      1                 2
#     /       / 1 + x'(y)  \
# min |  sqrt | ---------- |  dy
#     /       \    y       /
#     0
#

param n := 64;

param y {j in 0..n} := (j/n)^6;  # uneven spacing of the y partition is
				 # critical to get accurate solutions.
				 # Need lots of points at the beginning.

var x {j in 0..n};
var dxdy {j in 1..n} = (x[j] - x[j-1])/(y[j] - y[j-1]);
var f {j in 1..n} = sqrt( (dxdy[j]^2+1)/((y[j]+y[j-1])/2) );

minimize time: sum {j in 1..n} f[j]*(y[j]-y[j-1]) ;

subject to x0: x[0] = 0;
subject to xn: x[n] = 1;

#let {j in 0..n} x[j] := y[j]^(1/0.7);
let {j in 0..n} x[j] := y[j];

display time;

solve;

display time;

printf {j in 0..n}: "%10.5f %10.5f \n", x[j], -y[j] > bc.out;