Simple harmonic motion, RK4 method, matlab plot

Question

So I am having a issue plotting a simply harmonic motion of the form

$$\frac{d^2y}{dx^2}+\frac{k}{m}y=0$$

Using the RK4 method in matlab. The issue is that my code is not producing the expected plotted and I am not entirely sure, if it my RK4 that is wrong or my actual code that is wrong.

So the solution to the above equation is $y=3cos(\frac{k}{m}x)$, if i set k=1 and m=1, I should produce the following graph.

Using Demos

The graph that is being produced by my code is

here is the code which produces the given graph

clc;
clear all;
%Defing intial paramters
%dt=0.001;

h=0.001;
k=1;
m=1;
t=0:h:100;

x=zeros(length(t),1);
v=zeros(length(t),1);
x(1)=3;
v(1)=0;

%defing functions

f=@(t,x,v) v;
g=@(t,x,v) -k/m*x;





for n=1:(length(t)-1);
k1x=f(t(n),x(n),v(n));
k1v=g(t(n),x(n),v(n));

k2x=f(t(n)+0.5*h,x(n)+0.5*h*k1x,v(n)+0.5*h*k1v);
k2v=g(t(n)+0.5*h,x(n)+0.5*h*k1x+v(n)+0.5*h*k1v);

k3x=f(t(n)+0.5*h,x(n)+0.5*h*k2x,v(n)+0.5*h*k2v);
k3v=g(t(n)+0.5*h,x(n)+0.5*h*k2x,v(n)+0.5*h*k2v);

k4x=f(t(n)+h,x(n)+h*k3x,v(n)+h*k3v);
k4v=g(t(n)+h,x(n)+h*k3x+v(n)+h*k3v);

x(n+1)=x(n)+h/6*(k1x+2*k2x+2*k3x+k4x);
v(n+1)=v(n)+h/6*(k1v+2*k2v+2*k3v+k4v);


end

plot(t,x);

` I cannot understand to why my graph is being produced this way, I have gone through the code a couple of time, but cant see to find what is wrong. My only conclusion so far is that it either a variable I have missed or place wrong or my for stamen t is incorrect but cant see that being the case due to, using the equation for RK4 from a text book, which I have checked several times.

I have a very limited scoop when it comes to numerical programming and I am trying to get to better grips with i, and understand if there more efficient way of coding this but if possible I would like to understand what incorrect within my code, if possible. Any help would be much appreciated.

Answer 1

I do not know why that does not produce an error of insufficient arguments. You have two identical typos where instead of a comma there is a plus in the computation of k2v and k4v.

k2x=f(t(n)+0.5*h,x(n)+0.5*h*k1x,v(n)+0.5*h*k1v);
k2v=g(t(n)+0.5*h,x(n)+0.5*h*k1x+v(n)+0.5*h*k1v);

k4x=f(t(n)+h,x(n)+h*k3x,v(n)+h*k3v);
k4v=g(t(n)+h,x(n)+h*k3x+v(n)+h*k3v);