How to do Runge Kutta 4 with a second order ode?

I'm trying to do Runge Kutta with a second order ode, d2y/dx2+.5dy/dx+7y=0, with .5 step size from 0 to 5. My initial conditions are y'(0)=0 and y(0)=4.

This is my code so far:

s4=11; %Number of steps
h5=5/(s4-1); %Step size
x=linspace(0,5,s4);
y=[4; 0];

%Creating a for loop to define all the variables for Runge Kutta equation
for ii=1:length(x)-1
    k1(:,ii)=ode1(x(ii),y(:,ii));
    k2(:,ii)=ode1(x(ii)+.5*h5,y(:,ii)+k1(:,ii)*.5*h5);
    k3(:,ii)=ode1(x(ii)+.5*h5,y(:,ii)+k2(:,ii)*.5*h5);
    k4(:,ii)=ode1(x(ii)+h5,y(:,ii)+k3(:,ii)*h5);
    y(2,ii+1)=y(2,ii)+1/6*(k1(2,ii)+2*k2(2,ii)+2*k3(2,ii)+k4(2,ii))*h5;
    y(1,ii+1)=y(1,ii)+1/6*(k1(1,ii)+2*k2(1,ii)+2*k3(1,ii)+k4(1,ii));
end

I can't figure out why I'm not getting the right answer. Any ideas?

Answer 1

[James Tursa]

What is the signature for ode1? I can't seem to find it in the doc. Where is your derivative function implemented? I don't see it in your code. If you look at the following link for RK4, the k's represent function derivatives at various predicted points, not integrated values as you seem to be doing.

http://en.wikipedia.org/wiki/Runge_kutta

E.g., I would have expected to see a derivative function something like this in your code somewhere:

function yp = y_deriv(x,y)
  yp= [y(2);-7*y(1)-0.5*y(2)];
end

And this derivative function would be used to calculate the k's. But I don't see anything like this in your code. E.g.,

    k1(:,ii)=y_deriv(x(ii),y(:,ii));
    k2(:,ii)=y_deriv(x(ii)+.5*h5,y(:,ii)+k1(:,ii)*.5*h5);
    k3(:,ii)=y_deriv(x(ii)+.5*h5,y(:,ii)+k2(:,ii)*.5*h5);
    k4(:,ii)=y_deriv(x(ii)+h5,y(:,ii)+k3(:,ii)*h5);

Also looks like you are missing a *h5 on the last line. (you could combine the last two lines into one vector equation)

[Andrew]

The *h5 was actually the whole problem. Thanks a lot

Answer 2

[Diye]

try solving it using this method.

firstly split the second order ode into systems of first order odes. let, dy/dx=z then, dz/dx=-5z-7y

                 in Matlab 
function can %can is any name chosen for the script file

a0=[4 0]; [x, y]=ode45(@pol,[0:5:5],a0) y=a(:,1); z=a(:,2); plot(x,y) function dydx=pol(x,y) y=a(1); z=a(2); dydx=z; dydx=-5*z-7*y;