Solving ODE with singularity

Hey there,

I have trouble with solving an ODE. The code I'm running is pretty simple :

tspan = [0 10];
if l == 1
    y0 = [0 1];
elseif l == 0
    y0 = [1 0];
else
    y0 = [0 0];
end

[t,y] = ode45(@vortexODE, tspan, y0);

Where my vortexODE function is as follow :

function [ dydt ] = vortexODE( t,y )
    dydt = zeros(2,1);
    dydt(1) = y(2);
    dydt(2) = -(1/t)*y(2) + y(1)/(t^2) + y(1)^3 - t;
end

The thing is that when I try to solve the ODE, I get a vector of NaN, most probably due to the singularity at t=0 of my ODE. The thing is that the resulting function should be continuous around t=0 and is equal to 0 at t=0, and the initial conditions I have are only for t=0, and nothing more, and I have no idea how to work around this singularity to get ode45 actually provide me a result. Is there anything I could do to solve this ODE?

Answer 1

[Jan]

 the resulting function should be continuous around t=0 and is equal to 0 at t=0

No, this is not the case for y0=[0 1] and y0=[1 0]. And even if you fill the point at t=0 for y=[0,0], the integration fails due to the pole:

dydt = zeros(2,1);
dydt(1) = y(2);
if t == 0 && all(y==0)
    a = 0;
    b = 0;
else
    a = -(1/t)*y(2);
    b = y(1)/(t^2);
end
dydt(2) = a + b + y(1)^3 - t;

It does not look like the function can be integrated at t=0.