function [U_save,x,y,t_save]= ondas2D(c,t_f,M,tsteps);
%function [U_save,x,y,t_save]= ondas2D(c,t_f,M,tsteps);
%
%Resuelve utt - c2 (uxx +uyy) = 0 mediante un esquema explicito de segundo orden
%
%Datos de entrada:
% c: coeficiente c de la ecuacion
% M: numero de partes iguales en que se descompone cada intervalo [0,1]
% 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;
L=1;%Cuadrado [0,1]
h=L/M;
x=0:h:L;
y=0:h:L;
[X,Y]=meshgrid(x,y);
%%Condiciones iniciales
F=min(min(X,Y),min(1-X,1-Y));%La posicion inicial
G=zeros(size(X)); %La velocidad inicial
%Condiciones iniciales, sin los bordes:
F_red=F(2:M,2:M);
G_red=G(2:M,2:M);
%Versiones aplanadas
f=reshape(F_red,(M-1)^2,1);
g=reshape(G_red,(M-1)^2,1);
%%Datos de salida: en vez de guardar todos los datos intermedios,
%guardamos el vector u solo Nframes veces.
Nframes=30; marca=floor(tsteps/(Nframes-1));
%Utilizamos la matriz U para trabajar comodamente con el primer indice
%para la coordenada ’x’ y el segundo para la coordenada ’y’
U=zeros(M+1,M+1);
U_save=zeros(M+1,M+1,Nframes);
t_save=zeros(Nframes,1);
%Le ponemos las condiciones de frontera
%lado izquierdo
U_save(1,:,:)=zeros(1,M+1,Nframes);
%lado derecho
U_save(M+1,:,:)=zeros(1,M+1,Nframes);
%lado inferior
U_save(:,1,:)=zeros(M+1,1,Nframes);
%lado superior
U_save(:,M+1,:)=zeros(M+1,1,Nframes);
%guardamos la posicion de partida
U_save(:,:,1)=F;
t_save(1)=0;
%%El operador D
%Usamos el comando diag(vector,k) para crear una matriz tridiagonal
%Y el comando kron para poner los unos en los pares de indices
%que representan puntos vecinos del mallado
Vecinos = diag(-4*ones((M-1)^2,1))+...
kron(eye(M-1),diag(ones(M-2,1),1))+...
kron(eye(M-1),diag(ones(M-2,1),-1))+...
kron(diag(ones(M-2,1),1),eye(M-1))+...
kron(diag(ones(M-2,1),-1),eye(M-1));
D=(1/h^2)*Vecinos;
% Calculo del primer paso temporal
u_ant=f;
u = f + dt*g + 0.5*(dt^2)*(c^2)*D*f;
% 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(2:M,2:M)=reshape(u,M-1,M-1);
U_save(:,:,indice)=U;
t_save(indice)=tiempos(n);
end
end