# Maple derivation of the differential equation

Z:=sqrt(1+p^2)/(sqrt(2*g*y));
N:=diff(Z,y);
P:=diff(Z,p);
P:=subs({p=p(x),y=y(x)},P);
N:=subs({p=p(x),y=y(x)},N);
de:=factor(N-diff(P,x));
de:=subs({p(x)=diff(y(x),x),diff(p(x),x)=diff(y(x),x,x)},de);
desols:=dsolve({de,y(a)=A,y(b)=B},y(x));

s1:=desols[1];
eq1:=subs({y(x)=y},s1);
series((y-(1/2)*_C1)/sqrt(-y^2+y*_C1),y,3);

# see that in the limit when y goes to zero, arctan will have
# a negative argument going to -infinity, hence arctan=-Pi/2.
# Thus solve as a first condition 

_C1:=solve(1/2*_C1*Pi/2-_C2,_C1);

x:=1;
y:=1;
simplify(eq1);


# Matlab solution by finite differences

n=20;
h=1/(n+1);
e=ones(n,1);
D=spdiags([-e e],[-1 1],n,n)/2/h;
D2=spdiags([e -2*e e],[-1 0 1],n,n)/h^2;

a=1/2;
x=(h:h:1-h)';

y=a*x;
%y=a*exp(-x);
%y=zeros(size(x));
plot([0;x;1],-[0;y;a]);
bc=zeros(n,1);
bc(end)=a;

f=@(y) e+(D*y+bc/2/h).^2+2*(D2*y+bc/h^2).*y;
y=fsolve(f,y);
plot([0;x;1],-[0;y;a],'-o');

y=a*x;
w=0.002;
for i=1:1000
  r=e+(D*y+bc/2/h).^2+2*(D2*y+bc/h^2).*y;
  y=y+w*r;
  plot([0;x;1],-[0;y;a],'-o');
  pause
end;

x=20*(1+eps)
f=@(x) prod((1:50)-x)