I am trying to use Maple to apply Picard approximations to the following ODE:
$$ 2xyy' + x^2 - y^2 = 0,\ y(0) = 1 $$
I transformed it into the form $y' = f(x, y)$:
$$ y' = \frac {y^2 - x^2} {2xy} $$
Clearly, now the initial value problem is not in the domain of the function, but since I'm doing it numerically, I assumed it would be okay to take a «really small» value of $x$. Say, $x=10^{-10}$.
I followed the procedure described here (this is a copy-paste from my worksheet):
f := proc (x, y) options operator, arrow; (1/2)*(y^2-x^2)/(x*y) end proc;
x[0] := 10^(-10);
y[0] := 1
N := 3; # number of approximations
Now the tricky part.
for k to N do # in the code there actually was "k from 1 to N"
y[k] := unapply(y[0]+int(f(t, y[k-1](t)), t = x[0] .. x), x)
end do;
I do not understand why we have to use unapply here instead of a simple function definition like y[k] := x -> ... (Maple gives some weird errors with that), but okay. For now.
plot(y[1], 0 .. 1);
This works and produces a graph. However, for the second approximation
plot(y[2], 0 .. 1);
Maple just stops responding.
I tried the very same code for a different function
$$ f(x, y) = y^2 - x^2 $$
and it worked perfectly.
What am I doing wrong?
You are trying to use Maple to apply Picard approximations to the initial-value problem $$2xyy' + x^2 - y^2 = 0,\ y(0) = 1$$ This is just a wrong thing to do.
To begin with, there is no solution: the initial value is inconsistent with the ODE.
Second, you should consider what you are asking Maple to do. First, it is to plug $y_0=1$ into $\dfrac{y^2-x^2}{2xy} $ and integrate. Okay, that's a fair thing to ask: the result is $$y_1 = \frac{5}{4}-\frac{x^2}{4} +\frac12 \ln x$$ Then you asked it to plug $y_1$ into $\dfrac{y^2-x^2}{2xy} $ and integrate that. What were you thinking? Maple is not a magician to integrate such monsters.
When the right-hand side is a polynomial such as $y^2-x^2$, the iteration process produces polynomials, and everything works fine.
The form of your equation suggests that $u = y^2$ is a better function to look for: it satisfies the equation $xu'+x^2-u=0$ which is linear. Solutions are $u = Cx-x^2$ for arbitrary constant $C$. Hence, the original ODE has solutions $y = \sqrt{Cx-x^2}$. None of which satisfy the initial condition, of course.