Solving a system of 4 second order ODE's in matlab using ODE45

Question

I need to solve this system of second order equations using ODE45 in matlab I'm only familiar with using ODE45 with maybe one or two equations but not this many Here is what I have but I'm not sure how to correct it:

function second_oder_ode

t = 0:0.001:3;   % time scale
theta = pi/2;
phi = 0;

initial_t    = 0;
initial_dtds = 0;
initial_r    = 0;
initial_drds = 0;
initial_theta    = 0;
initial_dthetads = 0;
initial_phi    = 0;
initial_dphids = 0;

[t,x] = ode45( @(t,y) rhs(t,y,theta,phi), t, [initial_t initial_dtds initial_r initial_drds initial_theta initial_dthetads initial_phi initial_dphids] );

plot(t,x(:,1));
xlabel('t'); ylabel('r');
%%%%%%%%%%%%%% STATE VECTORS %%%%%%%%%%%%%%%%%%%%%
t = [t, dtds];
r = [r, drds];
theta = [theta, dthetads];
phi = [phi, dphids];
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function dxds=rhs(t,r, theta, phi)
    dxds_1 = t(2);
    dxds_2 = -((2/3)*R0^(2)*((3*H0/(2*cv))^(4/3))*y(1)^(1/3))*(y(4)^(2)+abs(y(3))^(2)*y(6)^(2)+abs(y(3))^(2)*(sin(y(5)))^(2)*y(8)^(2));
    dxds_3 = r(2);
    dxds_4 = -(4/(3*y(1)))*y(2)*y(4)+abs(y(3))*(y(6)^(2)+(sin(y(5)))^(2)*y(8)^(2));
    dxds_5 = theta(2);
    dxds_6 = -(2/abs(y(3)))*y(4)*y(6)-(4/(3*y(1)))*y(2)*y(6)+((sin(2*y(5)))*(y(8)^(2))/2);
    dxds_7 = phi(2);
    dxds_8 = -2*y(8)*((2/(3*y(1)))*y(2)+(cot(y(5)))*y(6)+(1/abs(y(3)))*y(4));

    dxdt=[dxdt_1; dxdt_2; dxdt_3; dxdt_4; dxdt_5; dxdt_6; dxdt_7; dxdt_8];

   end
end

Answer 1

This isn't a full solution, but is longer than a comment.

You seem to be using t theta and phi each for two purposes. For example, maybe change the timescale line to use s.

Looking at your odefun, here called rhs, I think it should have 2 arguments. The first argument should be s (even if it is not used) and the second argument should be a vector input where the values y(1:8) correspond to [t; t_s; r; r_s; theta; theta_s; phi; phi_s].

The expressions for dxds_# should only contain terms of y and no t, r, theta or phi.

You have typos on the line dxdt= which should contain dxds terms.

Below I've made some changes. I've not tested this; I don't have the constants R0, H0 or cv and I think it will be uniformly zero.

Also I've removed the first occurrences of theta and phi. I'm not sure how these are supposed to feature in rhs.

function second_oder_ode

s = 0:0.001:3;   % time scale
%theta = pi/2;
%phi = 0;

initial_t    = 0;
initial_dtds = 0;
initial_r    = 0;
initial_drds = 0;
initial_theta    = 0;
initial_dthetads = 0;
initial_phi    = 0;
initial_dphids = 0;

y0 = [initial_t, initial_dtds, initial_r, initial_drds, ...
      initial_theta, initial_dthetads, initial_phi, initial_dphids];

[s,x] = ode45( @(s,y) rhs(s,y), s, y0);

plot(s,x(:,1));
xlabel('s'); ylabel('t(s)');
%%%%%%%%%%%%%% STATE VECTORS %%%%%%%%%%%%%%%%%%%%%
%t = [t, dtds];
%r = [r, drds];
%theta = [theta, dthetads];
%phi = [phi, dphids];
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function dxds=rhs(s,y)
    dxds_1 = y(2);
    dxds_2 = -((2/3)*R0^(2)*((3*H0/(2*cv))^(4/3))*y(1)^(1/3))*(y(4)^(2)+abs(y(3))^(2)*y(6)^(2)+abs(y(3))^(2)*(sin(y(5)))^(2)*y(8)^(2));
    dxds_3 = y(4);
    dxds_4 = -(4/(3*y(1)))*y(2)*y(4)+abs(y(3))*(y(6)^(2)+(sin(y(5)))^(2)*y(8)^(2));
    dxds_5 = y(6);
    dxds_6 = -(2/abs(y(3)))*y(4)*y(6)-(4/(3*y(1)))*y(2)*y(6)+((sin(2*y(5)))*(y(8)^(2))/2);
    dxds_7 = y(8);
    dxds_8 = -2*y(8)*((2/(3*y(1)))*y(2)+(cot(y(5)))*y(6)+(1/abs(y(3)))*y(4));

    dxds=[dxds_1; dxds_2; dxds_3; dxds_4; dxds_5; dxds_6; dxds_7; dxds_8];

   end
end