Using Runge-Kutta algorithm to solve second order ODE

For this task I need to create a program which manually uses the RK algorithm to solve a second order ODE. I'm trying to get my head around it but I think I'm going at it in the wrong way.

My original function is with given values for m, F0 and w, so I have rearranged it to x'' = 50sin(8*pi*t)-x'-50x

I've then switched x' for y and rearranged to get 50sin(8*pi*t) = y' + y + 50x and this is what I'm using in the script. I need to plot x(t) and y(t) but at this point I'm just trying to find values for x and y over time.

The method I've used to get this far is below the code. Excuse my poor code. Thanks in advance.

PS I'm not allowed to use a function like ode45 to solve it for me.

close all; clear; clc;

t=0;x=0;y=0;
dy=myfunction(t,x,y)

x0=0;y0=0

h = 0.5;

for t=0:5.5
    
    dx1 = h*y
    dy1 = h*dy
    dx2 = h*(y+(dy1/2))
    dy2 = h*myfunction(t+(h/2),x+(dx1/2),y+(dy1/2))
    dx3 = h*(y+(dy2/2))
    dy3 = h*myfunction(t+(h/2),x+(dx2/2),y+(dy1/2))
    dx4 = h*(y+dy3)
    dy4 = h*myfunction(t+h,x+dx3,y+dy1)
    deltax = (dx1+(2*dx2)+(2*dx3)+dx4)/6
    deltay = (dy1+(2*dy2)+(2*dy3)+dy4)/6

    xth = x+deltax
    yth = y+deltay
    x=xth
    y=yth
    
    t = t+h
end
 

function [dy] = myfunction(t,x,y)
dy=(50*sin(8*pi*t))-(50*x)-y;
end

Method used (I swapped v for y):

Answer 1

[Torsten]

Can you take it from here ?

 %Solves y'' - (-exp(-B*t)-y+5*exp(-2*t)-2*exp(-(B+2)*t)+exp(-B*t)+t) = 0
 % t in [0 1]
 % y(0) = 1, y'(0) = -1
 % B = 4
 
  tstart = 0.0;
  tend = 1.0;
  h = (tend - tstart)/20;
  T = (tstart:h:tend).';
  Y0 = [1 -1]; 
  B = 4;
  f = @(t,y) [y(2) -exp(-B*t)-y(1)+5*exp(-2*t)-2*exp(-(B+2)*t)+exp(-B*t)+t];
  Y = runge_kutta_RK4(f,T,Y0);
  plot(T,Y)  
  
function Y = runge_kutta_RK4(f,T,Y0)
  N = numel(T);
  n = numel(Y0);
  Y = zeros(N,n);
  Y(1,:) = Y0;
  for i = 2:N
    t = T(i-1);
    y = Y(i-1,:);
    h  = T(i) - T(i-1);
    k0 = f(t,y);
    k1 = f(t+0.5*h,y+k0*0.5*h);
    k2 = f(t+0.5*h,y+k1*0.5*h);
    k3 = f(t+h,y+k2*h);
    Y(i,:) = y + h/6*(k0+2*k1+2*k2+k3);
  end
end              

[Daniel Jordan]

So can I change this line

f = @(t,y) [y(2) -exp(-B*t)-y(1)+5*exp(-2*t)-2*exp(-(B+2)*t)+exp(-B*t)+t];

into my own equation? Wouldn't it need to have t, x and y?

Also what do y(1) and y(2) mean in that line?

[Torsten]

For your function,

f = @(t,y) [y(2) F0/m*sin(omega*t)-c/m*y(2)-k/m*y(1)]

where

y(1) = y, y(2) = y'.

So you shouldn't work with x and y, but with a vector y=(y(1),y(2)) to make everything easier.

Of course, omega, F0, m and k have to be given values before.