function Bratzler_diffusion_eqn_students
% dimensionless parameters
sigma = 1;
alpha1 = 5;
alpha2 = 1;
C1s = 0.4;
C2s = 0.2;


t = 0; i = 0;
N = 50;
t_final = 3;
x = linspace(0,1,N); dx = x(2)-x(1);
dt = 0.5*dx^2;
C = zeros(1,N);


plot(x,C,'LineWidth',2);
tit = sprintf('time  = %f',t);
title(tit);
xlabel('x','FontSize',16); ylabel('LDL conc. C(x,t)','FontSize',16);
axis([0 1 0 1])

while t<t_final
    Cnew(2:N-1) = C(2:N-1) + dt*( (C(3:N)-2*C(2:N-1)+C(1:N-2))/dx^2 - sigma*(C(3:N)-C(1:N-2))/2/dx );
    Cnew(1) = ( Cnew(2)/dx + alpha1*C1s )/(sigma + 1/dx + alpha1 );
    Cnew(N) = ( -Cnew(N-1)/dx - alpha2*C2s )/( sigma - 1/dx - alpha2 );
    
    C = Cnew; t=t+dt; i=i+1;
    
    
    % if mod(i,30) == 0
    %     plot(x,C,'LineWidth',2);
    %     tit = sprintf('time  = %f',t);
    %     title(tit);
    %     xlabel('x','FontSize',16); ylabel('LDL conc. C(x,t)','FontSize',16);
    %     axis([0 1 0 1])
    %     pause(0.1)
    % end    
end
plot(x,C,'k+','LineWidth',2);
tit = sprintf('time  = %f',t);
title(tit);
xlabel('x','FontSize',16); 
ylabel('LDL conc. C(x,t)','FontSize',16);

