Error using ode45

I am trying to model the position of the Earth relative to the sun. I am using a vector form of the acceleration on a mass due to another mass. I have used this same formula in an Euler-Cromer scheme without error. However, when I try to implement it using the ode45 solver, I get this error:

 In an assignment  A(:) = B, the number of elements in A and B
must be the same.

 Error in orbit (line 10)
    dp(2) = -G * Msun * R / (sum(R.*R)^(3/2));

This is my function:

 function dp = orbit(t, P)
    %where p = [xi yi vxi vyi]   
     G = 6.673E-11/(1000^3);
    Psun = [0 0];
    Msun=1.988e30;
    R = Psun - P(2);
    M = 1.98892E30; 
    dp = zeros(4,1);
    dp(1)= P(2);
    dp(2) = -G * Msun * R / (sum(R.*R)^(3/2));

 end

And my script to call it:

 [T,Y] = ode45(@orbit,[0 12],[1.5e8 0 0 29.75]);

Answer 1

[Walter Roberson]

You define

Psun = [0 0];

so Psun is a vector. You then define

R = Psun - P(2);

so R is a vector.

Then when you have

dp(2) = -G * Msun * R / (sum(R.*R)^(3/2));

because you have that vector R on the right hand side, before the division, your result is going to be a vector. You are using /, the mrdivide operator, but (sum(R.*R)^(3/2)) is a scalar so R / (sum(R.*R)^(3/2)) is the same size as the vector R. So your overall expression is a vector, and you are trying to store that vector in the single location dp(2)

[Ibro Tutic]

Yea, I just realized that was causing the error. Using the vector form of that equation made using the Euler-Cromer method more straightforward for me. In this case, is there any work around besides doing everything in components?

Is there anyway to tell the program to expect a vector output and to store it accordingly?

[Walter Roberson]

dp(2:3) = -G * Msun * R / (sum(R.*R)^(3/2));

But then your ode would need an input vector of length 3 (because the output vector length must equal the input vector length)

You could put the two of them together by encoding as a single complex number, but then it would do complex integration on the term.