How to plot a simple wave equation using method of lines?

Hello, I am trying to build a simple wave equation solver that uses the method of lines and runge-kutta. The example I am testing is of the form u_tt=u_xx with initial conditions u(0,x)=sin(2pix), u_t(0,x) = 0. I know the exact solution to be cos(2pi*t)sin(2pi*x), whose range is from -1 to 1. However, running my solver gives me very large values, and I'm not sure where they are coming from. Here is my approximation (red), vs the exact solution (blue):

I think these values are way too big, even considering error. Can anyone point me to where I've gone wrong? I am not super familiar with the wave equation. Here is my program and Runge-Kutta solver:

n = 100; m = 100; tf = 1; 
x = linspace(0,1,n+1)'; xs = x(2:end-1); dx = 1/n; dt = 1/m;

f = @(x) sin((2*pi).*x); g = 0;                    % initial conditions
exact = @(t,x) (cos((2*pi).*t))*(sin((2*pi).*x));  % exact for comparison

% building second difference matrix
d  = kron(ones(n,1),[1,-2,1]); D2 = (1/dt^2) * spdiags(d, -1:1, n-1, n-1);

% build system of eqns
w0 = [f(xs);zeros(size(xs))];
t = linspace(0,tf,m+1)'; F = @(x,w) [w(n:end);D2*w(1:n-1)];
w = RK_solver(F,t,w0); w = reshape(w,n-1,(m+1)*2);

figure(1)
hold on
for i = 1:m+1
    %u = w(:,i);
    plot(x,exact(t(i),x),'b'); axis([0 1 -1 1])
    %pause(0.01);
end
hold off

figure(2)
hold on
for i = 1:m+1
    u = w(:,i);
    plot(x,[0;u;0],'r'); %hold on 
    %pause(0.01);
end

hold off

function w = RK_solver(F,t,c)

n = length(t)-1;               % number of subdivisions
m = length(c);                 % number of equations 
dt = t(2)-t(1);                % change in t 
w = zeros(m,n+1); 
w(:,1) = c;  % To store answers 

for i = 1:n
    k1 = dt.*F(t(i),          w(:,i));
    k2 = dt.*F(t(i)+(dt/2),   w(:,i)+(k1/2));
    k3 = dt.*F(t(i)+(dt/2),   w(:,i)+(k2/2));
    k4 = dt.*F(t(i)+dt,       w(:,i)+k3);
    k = (k1+(2*k2)+(2*k3)+k4)/6;
    w(:,i+1) = w(:,i) + k;
end

end

[Torsten]

Where do you plot the red solution curves in your code ? I can only see the exact solution in blue over time.

[Iris]

Sorry, I left that out on accident. It was just this loop:

for i = 1:m+1
    u = w(:,i);
    plot(x,[0;u;0],'r'); hold on 
    pause(0.01);
end

Answer 1

[Torsten]

n = 100; m = 100; tf = 1; dx = 1/n;
x = linspace(0,1,n+1)';
tspan = linspace(0,tf,m+1);
y0 = [sin(2*pi*x);zeros(n+1,1)];
%[t,y] = ode15s(@(t,y)fun(t,y,n,dx),tspan,y0);
y = RK_solver(@(t,y)fun(t,y,n,dx),tspan,y0);
y = y';
exact = cos(2*pi*tspan.').*sin(2*pi*x.');
figure(1)
hold on
for i = 1:m+1
plot(x,exact(i,:),'b')
end
hold off
axis([0 1 -1 1])

figure(2)
hold on
for i = 1:m+1
plot(x,y(i,1:n+1),'r')
end
hold off
axis([0 1 -1 1])

function dy = fun(t,y,n,dx)
  u = y(1:n+1);
  v = y(n+2:2*(n+1));
  dy = zeros(2*(n+1),1);
  dy(1:n+1) = v;
  dy(n+2:2*(n+1)) = [0;(u(3:n+1)-2*u(2:n)+u(1:n-1))/dx^2;0];
end

function w = RK_solver(F,t,c)
  n = length(t)-1;               % number of subdivisions
  m = length(c);                 % number of equations 
  dt = t(2)-t(1);                % change in t 
  w = zeros(m,n+1); 
  w(:,1) = c;  % To store answers 
  for i = 1:n
    k1 = dt.*F(t(i),          w(:,i));
    k2 = dt.*F(t(i)+(dt/2),   w(:,i)+(k1/2));
    k3 = dt.*F(t(i)+(dt/2),   w(:,i)+(k2/2));
    k4 = dt.*F(t(i)+dt,       w(:,i)+k3);
    k = (k1+(2*k2)+(2*k3)+k4)/6;
    w(:,i+1) = w(:,i) + k;
  end
end

[Iris]

Thank you. Do you know what was wrong with my Runge Kutta function, though? I am trying to not use a built in function.

[Torsten]

My guess is there is something wrong with the system of differential equations you pass to the Runge-Kutta function. The integrator looks correct (see above for a solution with your Runge-Kutta code, but my ODE system).

[Iris]

Ah okay. Thanks!