%Resumen: la idea del programa es crear una matriz de pares (tiempo,
%espacio) y en cada uno de esos puntos calcular la aproximación de la
%solución mediante un bucle for que realiza la suma de términos de la serie
%calculada

%creamos la malla en  la que introduciremos las evaluaciones 
evaluaciones2 = zeros(100,100); 
%quedan impuestas los condiciones frontera por ser una matriz de ceros ya
% que posteriormente no modificamos ni la primera ni la última columna

%limite temporal y vectores de tiempos y espacios
T = 2;
t = linspace(0,T,100); 
x = linspace(0,1,100); 

%funcion de la condicion inicial
cond =@(y) exp(-100.*(y-1/2).^2);
cond2 =@(y) (200.*y-100).*exp(-100.*(y-1/2).^2);


%en el siguiente bucle for, para cada k calculamos el coeficiente de
%Fourier y sumamos  el correspondiente término de la serie en cada punto de
% la malla. Para el cálculo de la integral que nos da el coeficiente de
% Fourier usamos la regla del trapecio con 1000 puntos en [0,1]
for k=1:50 
    %puntos
    N=1000;  
    %extremos del intervalo
    a=0; b=1;            
    h=(b-a)/N; 
    %partición del intervalo
    u=a:h:b;           
    %f es la funcion que queremos integrar
    f=(cond(u).*cos(k.*pi.*u))'; 
    g=(cond2(u).*cos(k.*pi.*u))';
    %vector de pesos
    w=ones(N+1,1);                 
    w(1)=1/2; w(N+1)=1/2;
    %m es el valor de la integral
    m=h*w'*f; 
    n=h*w'*g;
    %calculo final de los coeficientes del seno 
    coef = 2*m;   
    coef2 = 2*n;
    %sumamos los 100 primeros términos de la serie en cada punto de la
    % malla con los nuevos coeficientes
    for j =1:100
        for i = 1:100
                evaluaciones2(i,j) = evaluaciones2(i,j) + (coef.*cos(k*pi*t(i)) + coef2/(k*pi).*sin(k*pi*t(i))).*cos(k.*pi.*(x(j))); 
        end
    end
end

%plot en 3D
surf(t,x,evaluaciones2') 
xlabel('t')
ylabel('x')
zlabel('u(x,t)')