MATLAB: Runge-Kutta function with a second order ODE

differential equationsrunge kutta

Hello, I have this section of code I am trying to work with. I need it to be able to accept a matrix form of a second order derrivative. I'll post an example below. 
function [derriv_value] = FunctionC(x,y)
%Function that contains the derrivative value
derriv_value = [y(2); -9*y(1)+sin(x)]; % y(1) = y , y(2) = v
end
This is my function I am calling into my Runge-Kutta function. It is a second order ODE. I need my Runge-Kutta to be able to accept it, but I am not sure how. I tried altering how the inputs to the equation are formatted but nothing has worked. Here is the Runge-Kutta code.
function [x, yvecb] = MyVec_Function2(F,h,x0,x1,y0,y1)
% Note that F function expression is defined via Function Handle: F = @(x, y)(x+y)               
% change the function as you desire
                                             % step size (smaller step size gives more accurate solutions)
x = x0:h:x1;                                          % x space
yvecb(:,1)= y0;                         % Memory allocation
y(y0) = y1;                                       % initial condition
for i=1:(length(x)-1)                             
% i=1:(length(x)-1)                              % calculation loop
    k1 = F(x(i),y(:,i));
    k2 = F(x(i)+0.5*h,y(i)+0.5*h*k1);
    k3 = F((x(i)+0.5*h),(y(:,i)+0.5*h*k2));
    k4 = F((x(i)+h),(y(:,i)+k3*h));
    y(:,i+1) = y(:,i) + (1/6)*(k1+2*k2+2*k3+k4)*h;  % main equation
end
figure, plot(x, y)  % To see the solution results
end

Best Answer

I'm still confused about your initial conditions and what y0 and y1 are. But let's back up a bit.
For a 2nd order ODE, your state vector will have 2 elements. E.g., if y is your state then y(1) could be position and y(2) could be velocity. This looks like the setup you have in your FunctionC( ) above. So for that case, the states in your RK4 solver have to reflect this and be 2-element vectors, not scalars. So that y pre-allocation needs to have 2 rows, not 1. And all the downstream indexing needs to account for this as well. E.g., something like this if we assume that those y0 and y1 values are the initial values at the x0 time (CAUTION: Not checked):
function [x, y] = FunctionBeta_Executor(F,h,x0,x1,y0,y1)
% Note that F function expression is defined via Function Handle: F = @(x, y)(x+y)               
% change the function as you desire
                                             % step size (smaller step size gives more accurate solutions)
x = x0:h:x1;                                          % x space
y = zeros(2,length(x)); % CHANGED                            % Memory allocation
y(1,1) = y0;  % CHANGED                                       % initial condition
y(2,1) = y1;  % NEW
for i=1:(length(x)-1)                             
% i=1:(length(x)-1)                              % calculation loop
    k1 = F( x(i)         , y(:,i)            ); % CHANGED



    k2 = F( x(i) + 0.5*h , y(:,i) + 0.5*h*k1 ); % CHANGED
    k3 = F( x(i) + 0.5*h , y(:,i) + 0.5*h*k2 ); % CHANGED
    k4 = F( x(i) +     h , y(:,i) +     h*k3 ); % CHANGED
    y(:,i+1) = y(:,i) + (1/6)*( k1 + 2*k2 + 2*k3 + k4 )*h;  % main equation % CHANGED
end
figure, plot(x, y(1,:))  % To see the solution results % CHANGED
end
Note also how much easier it is to read (and debug) when you use spacing!
Related Question