NaN values appearing using ode45

I'm trying use to ode45 to solve a system of second order differential equations, but I keep getting NaN values. Here's the code:

syms xA(t) yA(t) xB(t) yB(t)

%---------------------------------------------------------------------

G = 1;

mA = 1;
mB = 1;

xA0 = -1;
dxAdt0 = 0;
yA0 = 2;
dyAdt0 = 0;
xB0 = 2;
dxBdt0 = 0;
yB0 = -1;
dyBdt0 = 0;

%---------------------------------------------------------------------

Initial = [xA0, dxAdt0, yA0, dyAdt0, xB0, dxBdt0, yB0, dyBdt0];

eq1 = diff(xA,2) == G*mB*(1/(((xB-xA)^2+(yB-yA)^2)^(3/2)))*(xB-xA);
eq2 = diff(yA,2) == G*mB*(1/(((xB-xA)^2+(yB-yA)^2)^(3/2)))*(yB-yA);
eq3 = diff(xB,2) == -G*mA*(1/(((xB-xA)^2+(yB-yA)^2)^(3/2)))*(xB-xA);
eq4 = diff(yB,2) == -G*mA*(1/(((xB-xA)^2+(yB-yA)^2)^(3/2)))*(yB-yA);

V = odeToVectorField([eq1,eq2,eq3,eq4]);
M = matlabFunction(V, 'vars', {'t', 'Y'});

[t,sol] = ode45(M, 0.1:1:5, Initial);

sol

sol = 5×8

-1 0 2 0 2 0 -1 0 NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN

Answer 1

[Walter Roberson]

Your initial conditions are incompatible with those equations -- with the divisions and those initial conditions you get divisions by 0.

It is tricky to find initial conditions that are valid, but the below shows an example that is able to integrate sometimes.

format long g
syms xA(t) yA(t) xB(t) yB(t)

%---------------------------------------------------------------------

G = 1;

mA = 1;
mB = 1;

xA0 = -1;
dxAdt0 = 0;
yA0 = 2;
dyAdt0 = 0;
xB0 = 2;
dxBdt0 = 0;
yB0 = -1;
dyBdt0 = 0;

%---------------------------------------------------------------------

Initial = [xA0, dxAdt0, yA0, dyAdt0, xB0, dxBdt0, yB0, dyBdt0];

eq1 = diff(xA,2) == G*mB*(1/(((xB-xA)^2+(yB-yA)^2)^(3/2)))*(xB-xA);
eq2 = diff(yA,2) == G*mB*(1/(((xB-xA)^2+(yB-yA)^2)^(3/2)))*(yB-yA);
eq3 = diff(xB,2) == -G*mA*(1/(((xB-xA)^2+(yB-yA)^2)^(3/2)))*(xB-xA);
eq4 = diff(yB,2) == -G*mA*(1/(((xB-xA)^2+(yB-yA)^2)^(3/2)))*(yB-yA);

[eqn, vars] = reduceDifferentialOrder([eq1, eq2, eq3, eq4], [xA, yA, xB, yB])

eqn = 

vars = 

[M, F] = massMatrixForm(eqn, vars);

f = M\F

f = 

odefun = odeFunction(f,vars)

odefun = function_handle with value:
    @(t,in2)[in2(5,:);in2(6,:);in2(7,:);in2(8,:);-(in2(1,:)-in2(3,:)).*1.0./(in2(1,:).*in2(3,:).*-2.0-in2(2,:).*in2(4,:).*2.0+in2(1,:).^2+in2(3,:).^2+in2(2,:).^2+in2(4,:).^2).^(3.0./2.0);-(in2(2,:)-in2(4,:)).*1.0./(in2(1,:).*in2(3,:).*-2.0-in2(2,:).*in2(4,:).*2.0+in2(1,:).^2+in2(3,:).^2+in2(2,:).^2+in2(4,:).^2).^(3.0./2.0);(in2(1,:)-in2(3,:)).*1.0./(in2(1,:).*in2(3,:).*-2.0-in2(2,:).*in2(4,:).*2.0+in2(1,:).^2+in2(3,:).^2+in2(2,:).^2+in2(4,:).^2).^(3.0./2.0);(in2(2,:)-in2(4,:)).*1.0./(in2(1,:).*in2(3,:).*-2.0-in2(2,:).*in2(4,:).*2.0+in2(1,:).^2+in2(3,:).^2+in2(2,:).^2+in2(4,:).^2).^(3.0./2.0)]

[t,sol] = ode45(odefun, 0.1:1:5, Initial);
sol

sol = 5×8

-1 0 2 0 2 0 -1 0 NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN

Ic = Initial+rand(size(Initial))/100

Ic = 1×8

-0.994089016105666 0.000472176043822199 2.00652394187924 0.00024701442291673 2.00571327100262 0.0026066757978841 -0.993910275708862 0.00333630442570037

[t,sol] = ode45(odefun, 0.1:1:5, Ic);

sol

sol = 5×8

-0.994089016105666 0.000472176043822199 2.00652394187924 0.00024701442291673 2.00571327100262 0.0026066757978841 -0.993910275708862 0.00333630442570037 14.1922600620998 25.4394580170961 -12.1680221410323 -25.4327958464059 102.902931015033 197.617243291868 -101.891128019738 -197.611300311645 117.095179934608 223.056679805197 -114.059139018246 -223.044074654283 102.902916943436 197.617216135336 -101.891113948142 -197.611273155112 219.998096294054 420.6738948135 -215.950252382397 -420.655346682364 102.902915977113 197.617214270444 -101.891112981818 -197.611271290221 322.901012074589 618.291108704573 -317.841365167638 -618.266617593213 102.902915628489 197.617213597641 -101.891112633194 -197.611270617418