"Grzegorz Knor" <gknor@o2.pl> wrote in message <i5js70$mgo$1@fred.mathworks.com>...
> Hi,
> I have the following code:
> 
> function test_ode
> del = pi;
> g = 0.6;
> T=1/(del*(1g));
> t = 0:.01:T;
>
> [tt yy] = ode45(@test_fun,t,0)
>
> function dy = test_fun(t,y)
> dy = del*(1y).^g;
> end
> end
> 
>
> Thi code solve the equation:
> dy/dt = del*(1y)^g
> Which analytical solution is:
> y(t) = 1((1g)(del*tc1))^(1/(1g))
> Having regard to the initial condition y(0) = 0:
> y(t) = 1(del*t*(g1)+1)^(1/(1g))
>
> Now supose that
> del = pi;
> g = 0.6;
> and time is:
> T=1/(del*(1g));
> t = 0:.01:T;
>
> For this vector t the analytical solution is real. So why the output from ode45 is complex?
>
> Best wishes
> Grzegorz
          
You are expecting too much of ode45. As y approaches one, as it must do at the endpoint T you have defined, it causes 1y to approach zero. Only a tiny error in the computed y at this point can cause 1y to go into the negative range and then you will get complex values for the next y. You had better stop a little before T if you don't want to run into that kind of trouble.
As you can see in your analytic solution, the point where y becomes one is a point of singularity and infinitely many branches in the complex plane can be followed from there. That should give you a clue that this region is to be avoided in what you are doing.
Roger Stafford
