Infinite while loop for Euler's method in order to solve ode

I am trying to solve an ode with Euler's method using while loop but i get into an infinite loop for my ymax value. If i change it to a lower value i don't have this problem. What is happening?

clc, clear all, close all

%--- Solving the diff. eqn. y'(t) = y(t) - c*y^2(t) with c = 0.5 using
%--- Euler's method in while loop

c = 0.5;
y0 = 0.1;
f = @(y) y - c*y.^2;
h = 0.1;
ymax = 1/c; % ymax is obtained by setting y' = 0 and solving for y, here is ymax = 2

t = 0;
y = 0.1;
[t,y]

while y < ymax % here is ymax = 2 and i get infinite loop. If, for example, i write y < 1 there is no problem, but i have to use ymax 2
    y = y + h*f(y);
    t = t + h;
    [t,y]
end

Answer 1

[John D'Errico]

I think you misunderstand what is happening. Your ODE is

y' = y - y^2/2

ymax is found when y' = 0, and so you find that ymax = 2. All good so far. But, WHEN does that happen? I'll solve the ode here:

syms y(t)
c = 0.5;
y0 = 1;
ODE = diff(y) == y - c*y^2;
yexact = dsolve(ODE,y(0) == 0.1)

yexact = 

fplot(yexact,[0,10])

Now, when does y(t) exceed 2? Never, of course. But does it EVER reach 2 EXACTLY? NO. Ony at t==infinity. And infinity is a long way away. But your test in the while triggers only at y==2 (or greater, which can never happen.)

limit(yexact,t,inf)

ans = 

2

while y < ymax % here is ymax = 2 and i get infinite loop. If, for example, i write y < 1 there is no problem, but i have to use ymax 2

And of course, if you change ymax to some lower value, it stops. Should that surprise you in the least bit? NO! Of course not.

[Left Terry]

OMG, you are so right and i feel so dumb right now. I changed it and, of course, it works. (I added some new lines trying to put the results in an array and plotting them).

clc, clear all, close all

%--- Solving the diff. eqn. y'(t) = y(t) - c*y^2(t) with c = 0.5 using
%--- Euler's method in while loop

c = 0.5;
f = @(y) y - c*y.^2;
h = 0.1;
ymax = 1/c;

t(1) = 0; y(1) = 0.1;

while y < 1.999999999999999
    y(end+1) = y(end) + h*f(y(end));
    t(end+1) = t(end) + h;
end

plot(t,y)

[Left Terry]

@John D'Errico Hello. When i run your code on my computer, the plot pops out but the axes limits are both from 0 to 1 instead being the ones we see below. What is wrong?

syms y(t)
c = 0.5;
y0 = 1;
ODE = diff(y) == y - c*y^2;
yexact = dsolve(ODE,y(0) == 0.1)

yexact = 

fplot(yexact,[0,10])

●

[Walter Roberson]

When I run John's code, the plot limits are 0 to 10.

[Left Terry]

@Walter Roberson That's not the case for me and i don't know why.

[Walter Roberson]

@Left Terry

You are encountering a bug specific to R2016a. The FunctionLine object returned by the plotting is using the wrong limits. The work-around is to use

xlim([0 10]);
ylim([0 2.2]);

https://www.mathworks.com/support/bugreports/details/1383817

[Left Terry]

@Walter Roberson I tried that, it didn't work. I restarded matlab and it didn't work either.

[Walter Roberson]

@Left Terry

I tested the code in R2016a. When I did not use xlim(), the output was restricted to [0 1] when I fplot with upper bound exceeding 1 . When I used xlim(), the plot worked properly.