bvp4c gives error "The derivative function ODEFUN should return a column vector of length 4."

I'm trying to use bvp4c and it's giving me the error above even though my function vector is a column vector of length 4. What am I missing here?

Primary script is:

global mu_Earth mu_Moon m0_til T_til m_dot_til;

m_Moon=7.349*10^22; %Mass of Moon, kg
m_Earth=5.9736*10^24; %Mass of Earth, kg
D=3.84402*10^5; %Mean Earth-Moon distance and reference length, km
w=2*pi/(27.322*24*3600);
w_til=w; %Average angular rate of Moon about Earth and reference time, rad/s
M=m_Moon+m_Earth; %Reference Mass, kg
G=(6.67428*10^-11)/(1000^3); %Universal gravitational constant, (km^3)/(kg*s^2)
T=5*12.1*10^-6; %Thrust, kg*km/s^2
m0=180-((24.5/w)*5*8*10^-7); %Initial spacecraft mass, kg
m_dot=5*8*10^-7; %Propellant mass flow rate, kg/s

m_til_Moon=m_Moon/M; %Non-dimensional mass of Moon
m_til_Earth=m_Earth/M; %Non-dimensional mass of Earth
t_til_max=0.4;
t_til_vec=0:.0001:t_til_max;
mu_Earth=((G*M)/((D^3)*(w_til^2)))*m_til_Earth; %Non-dimensional gravitational constant of Earth
mu_Moon=((G*M)/((D^3)*(w_til^2)))*m_til_Moon; %Non-dimensional gravitational constant of the Moon
T_til=T/(M*D*(w_til^2)); %Non-dimensional thrust
m0_til=m0/M; %Non-dimensional initial spacecraft mass
m_dot_til=m_dot/(M*w_til); %Non-dimensional mass flow rate

p=(3*pi/180)*sin(100*t_til_vec);

solinit=bvpinit(t_til_vec,'Bruh_AEM_594_guess',p);

tic
sol=bvp4c('Bruh_AEM_594_bvpfcn','Bruh_AEM_594_bcfcn',solinit);
toc

Equation function is:

function dcccdt_til=Bruh_AEM_594_bvpfcn(t_til,ccc,p)
% ccc(1)=r; ccc(2)=r_dot; ccc(3)=r*theta_dot; ccc(4)=theta
global mu_Earth mu_Moon m0_til T_til m_dot_til;

dcccdt_til=[ccc(2);
    ((ccc(3)^2)/ccc(1))+(((-(mu_Earth*cos(ccc(4)))/(ccc(1)^2))-((mu_Moon*((ccc(1)*cos(ccc(4)))+1))/((1+(2*ccc(1)*cos(ccc(4)))+(ccc(1)^2))^(3/2)))+((T_til*cos(ccc(4)+(pi/2)-p))/(m0_til-(m_dot_til*t_til)))+mu_Moon+(2*((ccc(2)*sin(ccc(4)))+(ccc(3)*cos(ccc(4)))))+(ccc(1)*cos(ccc(4))))*cos(ccc(4)))+(((-(mu_Earth*sin(ccc(4)))/(ccc(1)^2))-((mu_Moon*ccc(1)*sin(ccc(4)))/((1+(2*ccc(1)*cos(ccc(4)))+(ccc(1)^2))^(3/2)))+((T_til*sin(ccc(4)+(pi/2)-p))/(m0_til-(m_dot_til*t_til)))-(2*((ccc(2)*cos(ccc(4)))-(ccc(3)*sin(ccc(4)))))+(ccc(1)*sin(ccc(4))))*sin(ccc(4)));
    (-(ccc(2)*ccc(3))/ccc(1))-(((-(mu_Earth*cos(ccc(4)))/(ccc(1)^2))-((mu_Moon*((ccc(1)*cos(ccc(4)))+1))/((1+(2*ccc(1)*cos(ccc(4)))+(ccc(1)^2))^(3/2)))+((T_til*cos(ccc(4)+(pi/2)-p))/(m0_til-(m_dot_til*t_til)))+mu_Moon+(2*((ccc(2)*sin(ccc(4)))+(ccc(3)*cos(ccc(4)))))+(ccc(1)*cos(ccc(4))))*sin(ccc(4)))+(((-(mu_Earth*sin(ccc(4)))/(ccc(1)^2))-((mu_Moon*ccc(1)*sin(ccc(4)))/((1+(2*ccc(1)*cos(ccc(4)))+(ccc(1)^2))^(3/2)))+((T_til*sin(ccc(4)+(pi/2)-p))/(m0_til-(m_dot_til*t_til)))-(2*((ccc(2)*cos(ccc(4)))-(ccc(3)*sin(ccc(4)))))+(ccc(1)*sin(ccc(4))))*cos(ccc(4)));
    ccc(3)/ccc(1)];
end

Bounary condition function is:

function res=Bruh_AEM_594_bcfcn(ccca,cccb,p)
res=[ccca(1)-0.79 %Initial radius
    ccca(2)-0.2492 %Initial radial velocity
    ccca(3)-0.3865 %Initial tangential velocity
    ccca(4)-750.6497 %Initial angle
    cccb(1)-0.924671070669276 %Final radius
    cccb(2)-0.497482202882055 %Final radial velocity
    cccb(3)-0.2410837887939 %Final tangential velocity
    cccb(4)-750.797649549483]; %Final angle
end

And initial guess function is:

function g=Bruh_AEM_594_guess(t_til)
g=[(0.3815*(t_til^3))+(0.0401*(t_til^2))+(0.2592*t_til)+0.7898 %Radius
    (3.3006*(t_til^3))-(0.7603*(t_til^2))+(0.3895*t_til)+0.2465 %Radial Velocity
    (0.7827*(t_til^3))-(0.4514*(t_til^2))-(0.3115*t_til)+0.3856 %Tangential Velocity
    (0.0347*(t_til^3))-(0.3181*(t_til^2))+(0.4914*t_til)+750.65]; %Angle
end   

This was just a test that I know should converge with the solution giving value of p at every time step being 0.

Answer 1

[Stephan]

Your equations for dcccdt(2) and dcccdt(3) return vectors of size 4001x1, due to the length of p. This is the problem - there has to be a scalar result for both of them:

To fix this you should check your equations and correct them.

[Justin Bruh]

That makes sense.

I want it to evaluate p at the corresponding value of t_til for each time step. So p=p(27) at t_til=t_til(27), how do I code that into the equation function?

EDIT: Nevermind I fixed that by replacing p with p((t_til*10000)+1). But now it gives the error message:

"The boundary condition function BCFUN should return a column vector of length 4005)"

I only care about the 8 boundary conditions I gave it, what additional boundary conditions is it looking for?