function HeatEqn

k=1; L=1; t_end=1; N=10; s=0.52;

dx=L/N; dt=s*dx^2/k;

u_old=zeros(N+1,1);
%u_new=zeros(N+1,1);

%Initial conditions
x=zeros(N+1,1);
for ix=0:N
    j=ix+1; %Matlab runs from 1 to N+1, not 0 to N, so we need to shift it
    x_j=ix*dx;
    if x_j < 0.5*L
        u_old(j,1)=0;
    end
    if 0.5*L <= x_j && x_j < 0.75*L
        u_old(j,1)=(0.25/(0.25*L))*(x_j-0.5*L); %Line
    end
    if x_j > 0.75*L
        u_old(j,1)=-(0.25/(0.25*L))*(x_j-L);
    end
    
    x(j,1)=x_j; %Saving the grid
end

t=0;
u_new=u_old;
while t< t_end
    for ix=1:N-1
        j=ix+1;
        %PDE:
        u_new(j,1)=u_old(j,1)+s*(u_old(j+1,1)-2*u_old(j,1)+u_old(j-1,1));
    end
    
    %BC:
    u_new(1,1)=0; u_new(N+1,1)=0;
    
    u_old=u_new; %Update for the next step
    
    t=t+dt;
end

hold on
figure(1)
plot(x,u_new,'b -')
print('Solution','-depsc',figure(1))