So the theory is straightforward. We have:$$\frac{\partial^2U}{\partial t^2}=c^2 \frac{\partial^2U}{\partial x^2}$$discretizing it gives:$$\frac{U(i+1,j)- 2U(i,j) + U(i-1,j)}{(\Delta t)^2} = c^2 \frac{U(i,j+1)-2U(i,j) + U(i,j-1)}{(\Delta x)^2}$$where $U(i,j)$ is $U(t_i,x_j)$ with $t_i=i*\Delta t$ and $x_j=j*\Delta x$;Now, coding it up in Matlab I have:    function [u_new] = Explicit_W(dx,dt,n,initialize_pos,initialize_vel)%InputsX=1;   %X: length of the space domainJ=X/dx; %J: Total # of length intervalsf=str2func(initialize_pos); %initial conditions: u(0,x)g=str2func(initialize_vel); %initial conditions: du/dt(0,x)s2=(dt/dx)^  2;c=2;%Initializing the system, %u_ancient, u_old and u_new are U(i-1,:),U(i-1,:) and U(i+1,:) respectively.u_ancient=f(J);u_old=    f(J)+ dt*g(J)+ 1/2*c^2*Tridiag(J+1,s2,-2*s2,s2)*f(J);u_new=    Tridiag(J+1,c^2*s2,2-2*c^2*s2,c^2*s2)*u_old  - u_ancient;for i=1:n-2, % iterate n times, 2 of which are done already!    %Iterating the system using the Explicit formulation    temp=u_new;    u_new=Tridiag(J+1,c^2*s2,2-2*c^2*s2,c^2*s2)*u_old - u_ancient;    u_ancient=u_old;    u_old=temp;endendwhere Tridiag returns a matrix with entries on the tri-diagonals and its defined to be:function [ A ] = Tridiag( N,r1,r2,r3 )%This function creates a sparse tridiagonal matrix with fixed values%on the Diagonals.%   comes in handy when solving PDEs! i=[1:N 2:N   1:N-1]; %row indicesj=[1:N 1:N-1 2:N  ]; %column indicess=[(r2)*ones(1,N) r1*ones(1,N-1) r3*ones(1,N-1)]; %nonzero elementsA=sparse(i,j,s); % create the matrixendand for f, g I use Sin2 and Zero respectively, which are defined to be: function [ u0 ] = Sin2(L)%TENT Summary of this function goes here%   Initializes the array to Sin(2*pi*x)u0=zeros(floor(L+1),1);X=1;x=(0:X/L:1)';for j=1:(L+1),      u0=sin(2*pi*x);endendand function [ u0 ] = Zero( L)%Summary of this function goes here%   Initializes the array to a zero!u0=zeros(floor(L+1),1);endNow, solving it analytically, if use f=sin2 and g=zeros, we know the answer should be:$$u(x,t)=\frac{\sin(2\pi(x+ct))+\sin(2\pi(x-ct)}{2}$$so the expected range is $\{-1,1\}$. but when I execute Explicit_W, say:Explicit_W(0.05,0.0010,30,'Sin2','Zero')I get funny results, ranging from -124 to 124! Here is a picture of the results:http://i.stack.imgur.com/wON52.pngSo it has the form of the solution, but something is scaling it up! I have checked the sub-functions used in this on first order problems so I know there is nothing wrong with them. Whatever is causing the problem is in the main function. Any idea why it's acting the way it does?UPDATE:So it was brought to my attention that I had forgotten to implement the boundary conditions, here $u(0,t)=u(1,t)=0.$ Silly mistake!So the modified section of the code is:u_ancient=f(J);u_old=    f(J)+ dt*g(J)+ 1/2*c^2*Tridiag(J+1,s2,-2*s2,s2)*f(J);u_new=    Tridiag(J+1,c^2*s2,2-2*c^2*s2,c^2*s2)*u_old  - u_ancient;c1=0;c2=0;%boundary values%enforcing the boundary conditionsu_new(1)=c1; u_new(end)=c2;u_old(1)=c1; u_old(end)=c2;for i=1:n-2,    %Iterating the system using the Explicit formulation    temp=u_new;    u_new=Tridiag(J+1,c^2*s2,2-2*c^2*s2,c^2*s2)*u_old - u_ancient;    u_new(1)=c1;  %boundary conditions    u_new(end)=c2;%boundary conditions    u_ancient=u_old;    u_old=temp;endThat obviously fixes the issue with the boundary points, but I am still getting the wrong picture!! -> http://i.stack.imgur.com/lrVH7.pngIt's Working!I had fixed the issues that was pointed out and the code works fine now! it reads:    function [u] = Explicit_W(dx,dt,n,initialize_pos,initialize_vel)%Inputs%h: length of spatial intervalsX=1;   %X: length of the space domainJ=X/dx; %J: Total # of length intervalsf=str2func(initialize_pos);g=str2func(initialize_vel);s=(dt/dx);c=2;%velocityr=(s*c)^2;c1=0;c2=0;%boundary values at the two endsu=Zero(J+2);u(1)=c1;u(end)=c2;%output solution with its ends fixedb=Zero(J);b(1)=c1;b(end)=c2;  %vector b enforces the boundary conditions%u_new and u_old are solution vectors with their ends cut off.u_old=f(J);%Initializing u(0)u_new=f(J)+ dt*g(J)+ 1/2*r*(Tridiag(J-1,1,-2,1)*f(J)+b);for i=1:n-2,    %Iterating the system using the Explicit formulation    temp=u_new;     u_new=r*(Tridiag(J-1,1,-2,1)*u_new+b)+2*eye(J-1)*u_new - u_old;    u_old=temp;endu(2:end-1)=u_new;end
Implementing Explicit formulation of 1D wave equation in Matlab
pde;finite element;matlab;finite difference;linear solver
There may be several issues at work here, but one problem that stands out is your implementation of boundary conditions. You're working on the interval $[0,1]$, and partitioning it into $J = 1/dx$ subintervals. Including nodes on the boundary, this gives you $J+1$ total points, which we can call $x_1,\ldots, x_{J+1}$, with $x_1 = 0$ and $x_{J+1} = 1$. Because you are imposing Dirichlet boundary conditions, the value of $U$ at $x_1$ and $x_{J+1}$ is fixed throughout the simulation, and thus these points should not be included in your update matrix. Your update is of the form$$U(i+1,:) = (2I + c^2s_2A)U(i,:) - U(i-1,:)$$where $s_2 = (dt/dx)^2$ and$$A = \left(\begin{array}{ccccc} -2 & 1 & & & \\ 1 & -2 & 1 & \\  & 1 & -2 & 1\\ &&&\ddots&&\\&&&1&-2\end{array}\right)$$is the standard Laplacian matrix. But currently, you are applying this update to the vector $U(i,1:J+1)$, when you should be applying it to $U(i,2:J)$. Put differently, you should only be storing in u_ancient, u_old, and u_new the value of $U$ at $x_2,\ldots,x_{J}$. In fact, this formulation of the problem -- and homogenous Dirichlet boundary conditions -- are encoded directly into the matrix $A$ that you are (implicitly) using. Observe that\begin{align*}AU(i,2:J) &= \left(\begin{array}{c} -2U(i,1) + U(i,2) \\ U(i,1) - 2U(i,2) + U(i,3) \\ \vdots\\U(i,J-1) - 2U(i,J)\end{array}\right)\\&=\left(\begin{array}{c} U(i,0) -2U(i,1) + U(i,2) \\ U(i,1) - 2U(i,2) + U(i,3) \\ \vdots\\U(i,J-1) - 2U(i,J) + U(i,J+1)\end{array}\right)\end{align*}Here the first step is just matrix multiplication, but the second step follows precisely because you have assumed homogenous Dirichlet boundary conditions to begin with, i.e. because $U(i,0) = U(i,J+1) = 0$. This hopefully explains why you needn't include the points $x_1$ and $x_{J+1}$ in your system. EDIT: Upon re-reading the question, it occurs to me that the issue outlined above shouldn't have altered the solution on account of the problem's explicit formulation, so I'm afraid this doesn't explain your aberrant results. ANOTHER EDIT: I think I found your problem: the way you've named your variables has led you to implement the time-stepping incorrectly.Note that after the initial assignment of u_ancient, u_old, and u_new, it is the variable u_new that holds the current value of $U$. In other words u_old is $U(i-1,:)$ and u_new is $U(i,:)$. So in the next step, you want to base the new value, $U(i+1,:)$, off of the current value, $U(i,:)$, and the value that is one step behind, $U(i-1,:)$. You never actually need the value that is two timesteps behind, which you have named u_ancient. In contrast, you've been performing the update$$U(i+1,:) = (2I+c^2s_2A)U(i-1,:) - U(i-2,:)$$and this is inconsistent (in the techinical sense) with the PDE.The fix is very simple. You don't need the variable u_ancient at all, for one. And your loop should be replaced byfor i=1:n-2    temp=u_new    u_new=Tridiag(J+1,c^2*s2,2-2*c^2*s2,c^2*s2)*u_new - u_old    u_old=temp;end