Solving a complicated system of ODEs

I'm trying to numerically solve a set of differential equations by using ode45 (I don't know if this is the most appropriate one). My code is below. The annoying thing is that this seems to run but gives invalid solutions. Also, the k parameter should be a variable, not a parameter with a specific value. I don't know how to implement this. It is certainly possible that my initial conditions are incorrect, however, this shouldn't really have an impact if the code runs or not I think.

I'm not fluent at all with this syntax, so bare with me for not understanding everything.

EDIT: I've added the problem I am solving below. Note that the transfer function depends on k and can be obtained through δ I believe. I'll have to check this.

clear; clc;


sol = ode45(@ode_fun, [0 1], [1; 0; 0; 0; 0; 1]);


function dydt = ode_fun(t, y)
    rhodm = 0.8;
    rho = 1.0;
    G = 0.000000000066743;
    H0 = 1; %Model dependent?
    Omr = 0.0000463501;
    Omm = 0.2514;
    Omk = 0;
    OmLa = 0.6847;
    k=1.0;

    theta0 = y(1);
    theta1 = y(2);
    Phi = y(3);
    delta = y(4);
    v = y(5);
    a = y(6);
    
    dydt = zeros(size(y));
    dydt(1) = -dydt(5) - k.*theta1;
    dydt(2) = -k./3.*(Phi-theta0);
    dydt(3) = -3.*dydt(5) - 1i.*k.*v;
    dydt(4) = 1i.*k.*Phi - dydt(6)/a.*v;
    dydt(5) = (4.*pi.*G.*a.^2.*(rhodm.*delta + 4.*rho.*theta0) - k.^2.*Phi ).*a/(3.*dydt(6)) - dydt(6)/a.*Phi;
    dydt(6) = H0 * a .* sqrt(Omr ./ a.^4 + Omm ./ a.^3 + Omk ./ a.^2 + OmLa);
end

[James Tursa]

Please post an image/pic of the differential equations you are solving. If k is supposed to be a variable, what are the differential equations governing it? Why isn't it one of your y elements?

[bozo]

@James Tursa I've updated with the differential equations

[Torsten]

@bozo

You forgot the 5th equation.

Answer 1

[Torsten]

You use dydt(5) and dydt(6) before they have been assigned their values:

dydt(1) = -dydt(5) - k.*theta1;

dydt(2) = -k./3.*(Phi-theta0);

dydt(3) = -3.*dydt(5) - 1i.*k.*v;

dydt(4) = 1i.*k.*Phi - dydt(6)/a.*v;

dydt(5) = (4.*pi.*G.*a.^2.*(rhodm.*delta + 4.*rho.*theta0) - k.^2.*Phi ).*a/(3.*dydt(6)) - dydt(6)/a.*Phi;

Use instead a modified order in the calculation of the derivatives:

    dydt(6) = H0 * a .* sqrt(Omr ./ a.^4 + Omm ./ a.^3 + Omk ./ a.^2 + OmLa);
    dydt(5) = (4.*pi.*G.*a.^2.*(rhodm.*delta + 4.*rho.*theta0) - k.^2.*Phi ).*a/(3.*dydt(6)) - dydt(6)/a.*Phi;
    dydt(1) = -dydt(5) - k.*theta1;
    dydt(2) = -k./3.*(Phi-theta0);
    dydt(3) = -3.*dydt(5) - 1i.*k.*v;
    dydt(4) = 1i.*k.*Phi - dydt(6)/a.*v;