clear;

l = 1;

function dy = boat(x,y)
  dy = (y - b * sqrt(x ^ 2 + y ^ 2)) / x
endfunction

function y = boat_an(x)
  y = l / 2 * ((x / l) ^ (1 - b) - (x / l ) ^ (1 + b))
endfunction

y0 = 0;
X = linspace(l,0.01,40);

b = 0.1;
plot(X,boat_an(X),'-r');
Y = ode(y0, l, X, boat);
plot(X,Y,'or');

b = 0.5;
plot(X,boat_an(X),'-g');
Y = ode(y0, l, X, boat);
plot(X,Y,'og');

b = 0.8;
plot(X,boat_an(X),'-b');
Y = ode(y0, l, X, boat);
plot(X,Y,'ob');

b = 1.0;
plot(X,boat_an(X),'-m');
Y = ode(y0, l, X, boat);
plot(X,Y,'om');

b = 1.5;
plot(X,boat_an(X),'-c');
Y = ode(y0, l, X, boat);
plot(X,Y,'oc');

b = 2.0;
plot(X,boat_an(X),'-y');
Y = ode(y0, l, X, boat);
plot(X,Y,'oy');

zoom_rect([0,0,1,2]);
xlabel("x");
ylabel("y");