function [u_save,x,t_save]=calor2(k, M, t_f, tsteps)
%function [u_save,x,t_save]=calor2(k, M, t_f, tsteps)
%
%resuelve la ecuacion del calor u_t = k u_xx + f en [0,1] con dato inicial 0
%y condiciones de contorno tipo Neumann homogeneas
%
%Datos de entrada:
% k: coeficiente k de la ecuacion
% M: numero de partes iguales en que se descompone [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;
h=1/M;
x=0:h:1;

%Datos iniciales
u=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,1),1)+diag(ones(M,1),-1);
Tridiag(1,2)=2;
Tridiag(M+1,M)=2;
D=(1/h^2)*Tridiag;

%El vector f
mitad=floor((M+1)/2);
f=[ones(mitad,1); zeros(M+1-mitad,1)];

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

%Bucle principal
I=eye(M+1);
A=(I+dt*k*D);
b=dt*f;
for n=1:tsteps
    u=A*u+b;
    %Guardamos los valores de u para algunos tiempos
    if mod(n,marca)==0
      indice=1+n/marca;
      u_save(:,indice)=u;
      t_save(indice)=tiempos(n);
    end
end
