title: ode45 Output Vector and Initial Condition Vector
body: <p>Here is the error message I am receiving upon running my code:</p>
  
  <pre><code>Error using odearguments (line 92)
  @(T,Y,P)THREE_BODY_GRAV(T,Y,[4*PI^2,1,3.0135E-6,3.693E-8]) returns a vector of length
  12, but the length of initial conditions vector is 8. The vector returned by
  @(T,Y,P)THREE_BODY_GRAV(T,Y,[4*PI^2,1,3.0135E-6,3.693E-8]) and the initial conditions
  vector must have the same number of elements.
  </code></pre>
  
  <p>I don't understand why the output vector would be of length 12. The problem I am trying to solve is the three body gravitational system, in which the earth and moon orbit about the stationary sun. Here is my preliminary work (defining the functions, etc):</p>
  
  <pre><code>syms x_SE y_SE x_SM y_SM x_EM y_EM v_Ex v_Ey v_Mx v_My Ms Me Mm G t
  
  r_SE = [x_SE;y_SE]
  r_SM = [x_SM;y_SM]
  x_EM = x_SE - x_SM
  y_EM = y_SE - y_SM
  r_EM = [x_EM;y_EM]
  
  eqs_1 = -(Ms*G*x_SE)./(r_SE).^3-(Mm*G*x_EM)./(r_EM).^3
  eqs_2 = -(Ms*G*y_SE)./(r_SE).^3-(Mm*G*y_EM)./(r_EM).^3
  eqs_3 = -(Ms*G*x_SM)./(r_SM).^3-(Me*x_EM)./(r_EM).^3
  eqs_4 = -(Ms*G*y_SM)./(r_SM).^3-(Me*y_EM)./(r_EM).^3
  
  eqs = [eqs_1;v_Ex;eqs_2;v_Ey;eqs_3;v_Mx;eqs_4;v_My]
  
  matlabFunction(eqs,'vars',{t,[v_Ex;x_SE;v_Ey;y_SE;v_Mx;x_SM;v_My;y_SM],...
                     [G,Ms,Me,Mm]},'file','three_body_grav')
  </code></pre>
  
  <p>Looking back at this, I now realize that I poorly named the variables; but let me attempt at describing each one. <code>r_SE</code> is the vector that points from the sun to the earth; <code>r_SM</code> is the vector that points between the sun and moon; <code>r_EM</code> is the vector that points from the earth to the moon; <code>eqs1</code> to <code>eqs4</code> are the accelerations of the earth and moon in the x and y directions; and <code>eqs</code> is just a column vector to store these unknown function. <code>v_Ex</code> and <code>v_Ey</code> are the x and y direction velocities of the earth.</p>
  
  <p>I then used <code>matlabFunction</code> to define, which it seemed to do properly. However, once I reached the point in which I wanted to approximate the unknown functions, so that I can get the positions as a function of time, I received the error I alluded to above. Here is what I ran:</p>
  
  <pre><code>[tout,yout] = ode45(@(t,Y,P)three_body_grav(t,Y,[4*pi^2,1,3.0135e-6,3.693e-8]),...
                      [0,1],[0;1;1;0;0;1;1;0],myoptions);
  </code></pre>
  
  <p>To my understanding, there should only be eight initial conditions: one for <code>v_Ex</code> (initial x velocity of earth), <code>x_SE</code> (initial positon of the earth), <code>v_Ey</code> (initial y velocity of the earth); <code>y_SE</code> (initial y position of the earth), <code>v_Mx</code> (initial x velocity of the moon), <code>x_SM</code> (initial x position of the moon), <code>v_My</code> (initial y velocity of the moon), <code>y_SM</code> (initial y position of the moon).</p>
  
  <p>Could some possible interpret my error?</p>
  
  <p>EDIT: I was able to fix it. Here is the updated code:</p>
  
  <pre><code>syms x_E y_E x_M y_M x_EM y_EM v_Ex v_Ey v_Mx v_My Ms Me Mm G t
  
  r_E = [x_E;y_E]
  
  mag_r_E = sqrt(r_E.'*r_E)
  
  r_M = [x_M;y_M]
  
  mag_r_M = sqrt(r_M.'*r_M)
  
  r_EM = r_E - r_M
  
  mag_r_EM = sqrt(r_EM.'*r_EM)
  
  eqs_1 = -(Ms*G*x_E)./(mag_r_E).^3-(Mm*G*(x_E-x_M))./(mag_r_EM).^3
  
  eqs_2 = -(Ms*G*y_E)./(mag_r_E).^3-(Mm*G*(y_E-y_M))./(mag_r_EM).^3
  
  eqs_3 = -(Ms*G*x_M)./(mag_r_M).^3-(Me*G*(x_E-x_M))./(mag_r_EM).^3
  
  eqs_4 = -(Ms*G*y_M)./(mag_r_M).^3-(Me*G*(y_E-y_M))./(mag_r_EM).^3
  
  eqs = [eqs_1;v_Ex;eqs_2;v_Ey;eqs_3;v_Mx;eqs_4;v_My]
  
  matlabFunction(eqs,'vars',{t,[v_Ex;x_E;v_Ey;y_E;v_Mx;x_M;v_My;y_M],[G,Ms,Me,Mm]},'file','three_body_grav')
  </code></pre>
  
  <p>The problem was that I was dividing by a vector; and in Newton's law of gravity, we are supposed to divide by the magnitude of the distance vector.</p>
  
  <p>However, now I am having another difficulty. When I use </p>
  
  <pre><code>[tout,yout] = ode45(@(t,Y,P) three_body_grav(t,Y,[4*pi^2,1,3.0135e-6,3.693e-8]),[0,1],[0;1;1;0;0;1;1;0],myoptions);
  </code></pre>
  
  <p>to solve the unknown functions, all but the first element in the output vector is a 'NaN.' For instance,</p>
  
  <pre><code> yout(:,1)
  
  ans =
  
       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
     NaN
     NaN
     NaN
     NaN
     NaN
     NaN
     NaN
     NaN
  </code></pre>