Second Order Diff Eq with ode45 in Matlab

Question

So I need to solve x''(t) = -x(t)^p with initial conditions x(0)= 0 and v(0) = x'(0) = v_o = 1. The value of the parameter p is 1.

This is what I have:

function [t, velocity, x] = ode_oscilation(p)

y=[0;0;0];
    % transform system to the canonical form

    function y = oscilation_equation(x,p)
        y=zeros(2,1);
        y(1)=y(2);
        y(2)=-(x)^p;
        %  to make matlab happy we need to return a column vector
        % so we transpose (note the dot in .')
        y=y.'; 
    end

    tspan=[0, 30]; % time interval of interest

    [t,velocity,x] = ode45(@oscilation_equation, tspan, 1); 

    t = y(:,1);
    xposition=y(:,3);
    velocity=y(:,2); 

end 

and this is the error message I receive:

ode_oscillation(1) Error using odearguments (line 91) ODE_OSCILLATION/OSCILATION_EQUATION must return a column vector.

Error in ode45 (line 114) [neq, tspan, ntspan, next, t0, tfinal, tdir, y0, f0, odeArgs, odeFcn, ...

Error in ode_oscillation (line 17) [t,velocity,x] = ode45(@oscilation_equation, tspan,1);

Answer 1

There's a few things going wrong here. First, from help ode45:

ode45 Solve non-stiff differential equations, medium order method.

[TOUT,YOUT] = ode45(ODEFUN,TSPAN,Y0) with TSPAN = [T0 TFINAL] integrates 
the system of differential equations y' = f(t,y) from time T0 to TFINAL 
with initial conditions Y0. 

Note that ode45 expects a function f(t,y), where size(t) == [1 1] for time and size(y) == [1 N] or [N 1] for solution values. Your oscilation_equation has the order of input arguments inverted, and you input a constant parameter p instead of time t.

Also, the initial conditions Y0 should have the same size as y; so size(y0) == [N 1] or [1 N]. You just have 1, which is clearly causing errors.

Also, your output arguments t, xposition and velocity will be completely ignored and erroneous, since y is not set as output argument from ode45, and most of all, their names do not correspond to ode_oscilation's output arguments. Also, their order of extracting from columns of y is incorrect.

So, in summary, change everything to this:

function [t, v, x] = ode_oscilation(p)

    % initial values
    y0 = [0 1];           

    % time interval of interest
    tspan =[0 30]; 

    % solve system 
    [t,y] = ode45(@(t,y) [y(2); -y(1)^p], tspan, y0);

    % and return values of interest
    x = y(:,1);
    v = y(:,2);

end