Solve (numerically) a second-order ODE can someone help me in solving the following ODE:
$\sqrt y = ax - bx^2 (1-x)^2 + c x^2 (1-x)^2 y''$
$y(0)=0, y(1)=a^2$
The main problem is that boundary conditions directly follow from the equation itself. Also if you know how to correct it with a proper substitution of variables would be very appreciated.
This is an equation that emerges in an exceptionally important question on the verge of Economics and Mathematics (Multi-Armed Bandits).
 A: So I wrote some MATLAB code (You can also run it with Octave). I would prefer to use Chebfun but it sometimes causes problems with Octave and you might not have MATLAB. Here it goes

h = 0.001; x = 0:h:1;
a=2; b=2; c=1;
y = zeros(size(x));
y(1) = 0; y(end) = a^2;

objfun = @(x,y) (a*x(2:end-1) - ...
                b*x(2:end-1).^2.*(1-x(2:end-1)).^2 + ...
                c*x(2:end-1).^2.*(1-x(2:end-1)).^2.*(y(1:end-2) - 2*y(2:end-1) + y(3:end))./h^2)).^2;
ynew = y;
for i=1:50
  ynew(2:end-1) = objfun(x,y);
  norm(y-ynew)
  y = ynew;
endfor

plot(x,y)

Explanation:
First, I am approximating $y^{\prime\prime}(x) = (y(x-h)-2y(x)+y(x+h))/h^2$.
Second, I am considering the problem as a fixed point problem, i.e. find $y(x)$ such that $$\sqrt{y(x)} = ax - bx^2 (1-x)^2 + c x^2 (1-x)^2 (y(x-h)-2y(x)+y(x+h))/h^2.$$ Further, I am writing it as $$y(x) = (ax - bx^2 (1-x)^2 + c x^2 (1-x)^2 (y(x-h)-2y(x)+y(x+h))/h^2)^2.$$
Then I am using Picard iterations (that is the for loop). I am iterating a fixed number of steps but the difference between consecutive iterates is going to zero.
HOWEVER, I haven't analyzed the method! It may not be convergent, so be careful and check the result to see if it makes sense. I also haven't tested it thoroughly. Let me know if you have further questions, or if you encounter any problems.
