1D wave equation using Finite difference method MATLAB

Question

I have the wave equation $$u_{tt} = 4 u_{xx}$$ with the boundary conditions $$u(0,t) = u(L,t) = 0\,,\quad x \leq 0 \leq 2\pi \,,\quad t\geq 0$$ and initial conditions $$\begin{align} &u(x,0)=\begin{cases}1 & \text{for} \quad \pi-1 \leq x \leq \pi +1 \\0 &\text{otherwise} \end{cases} \\[0.5em] &u_{t}(x,0) = x^2sin^2(4x) \end{align}$$

I'm trying to implement this problem on MATLAB by the finite difference method and by using the surf function to plot it as a 3D wave; however, the problem I'm having is how to code the first initial condition.

The function I am using is

function [u, q] = Wave(f1,f2,g0,g1,xspan,tspan,nx,nt,alpha)
x0 = xspan(1)
xf = xspan(2)
t0 = tspan(1)
tf = tspan(2)
dx = (xf - x0)/nx;
dt = (tf-t0)/nt;
x = [0:nx]'*dx;
t = [0:nt]*dt;
q = alpha*(dt/dx)^2;
q1 = q/2;
q2 = 2*(1-q);
u(:,1) = f1(x);
for k = 1:nt+1, u([1 nx+1],k) = [g0(t(k)); g1(t(k))]; end
u(2:nx,2) = q1*u(1:nx-1,1) + (1-q)*u(2:nx,1) + q1*u(3:nx+1,1) + dt*f2(x(2:nx));
for k = 3:nt+1
    u(2:nx,k) = q*u(1:nx-1,k-1) + q2*u(2:nx,k-1) + q*u(3:nx+1,k-1) - u(2:nx,k-2);
end

surf(t,x,u) 
xlabel('t')
ylabel('x')
zlabel('u(x,t)')
end

and I know that within the function, the variable f1 is the one controlled by the first initial condition. In the function, I don't know how to incorporate the initial condition into the method. I assume that a for loop and an if statement are to be used but anything I tried doesn't work.

The script that I am using to plot the graph is

clc
clear all
f1 = @(x) ////;
f2 = @(x) (x.^2).*((sin(4.*x)).^2);
g0 = @(t) 0;
g1 = @(t) 0;
xspan = [0,2*pi];
tspan = [0,1];
nx = 20;
nt = 40;
alpha = 4;
[u,r] = Wave(f1,f2,g0,g1,xspan,tspan,nx,nt,alpha);

As seen on the script where f1 = @(x) ////;. I'm unsure of what to put here for the initial condition. Any help would be very much appreciated.

Answer 1

Check the interval of x and then covert logical to double.

f1 = @(x) double(pi-1 <= x && x <= pi+1);

Or if the initial condition is not 1 and 0, write a function

function u = f1(x)
    u = zeros(size(x));
    is_one = pi-1 <= x && x <= pi+1;
    u(is_one) = 1;
    % u(~is_one) = 0;  not required, just for your reference
end

Edit:

$u(x, 0)$ is set by

u(1,:) = f1(x);

Then you used it to calculate $u(x,\Delta t)$, i.e. u(:,2) in this line

u(2:nx,2) = q1*u(1:nx-1,1) + (1-q)*u(2:nx,1) + q1*u(3:nx+1,1)...
            + dt*f2(x(2:nx));

You have done it, haven't you?