function [u_save,x,t_save]= ondas(c, L, t_f, M,tsteps);
%function [u_save,x,t_save]= ondas(c, L, t_f, M,tsteps);
%
%Resuelve u_tt - c2 u_xx = 0 mediante un esquema explicito de segundo orden
%
%Datos de entrada:
% c: coeficiente c de la ecuacion
% L: longitud del intervalo
% M: numero de partes iguales en que se descompone [0,L]
% t_f: instante de tiempo hasta el que calculamos la solucion
% tsteps:numero de partes iguales en que se descompone [0,t_f]

%%Mallado
dt=t_f/tsteps;
tiempos=0:dt:t_f;
h=L/M;
x=0:h:L;

%%Condiciones iniciales
f=(sin(pi*x)+0.5*sin(3*pi*x))’;
g=zeros(M+1,1);

%%El operador D
%Usamos el comando diag(vector,k) para crear una matriz tridiagonal
Tridiag = diag(-2*ones(M-1,1)) + diag(ones(M-2,1),1) + diag(ones(M-2,1),-1);
D=(1/h^2)*Tridiag;

%%Datos de salida: en vez de guardar todos los datos intermedios,
%guardamos el vector u solo Nframes veces.
Nframes=11;
marca=floor(tsteps/(Nframes-1));
u_save=zeros(M+1,Nframes);
t_save=zeros(1,Nframes);
%Le ponemos las condiciones de frontera
u_save(1,:)=zeros(1,Nframes);
u_save(M+1,:)=zeros(1,Nframes);
%guardamos la posicion de partida
u_save(:,1)=f;
t_save(1)=0;

% Posicion de partida
u=f(2:M);
% Calculo del primer paso temporal
u_ant=u;
g_red=g(2:M);
u = u + dt*g_red + 0.5*(dt^2)*(c^2)*D*u;

% Iteracion para el resto de pasos temporales
%Bucle principal
for n=1:tsteps
    temp=u;
    u = 2*u - u_ant + (dt^2)*(c^2)*D*u;
    u_ant=temp;
    %Guardamos los valores de u para algunos tiempos
    if mod(n,marca)==0
        indice = 1 + n/marca;
        u_save(2:M,indice)=u;
        t_save(indice)=tiempos(n);
    end
end
