function fig229

%  plot of solution of 
%  u_t + u u_x = ep*u_xx  for  x0 < x < x1  
%  u(x,0) = step(u_1,u_2) where u_1 > u_2

clf

% get(gcf)
%set(gcf,'Position', [818 364 560 420]);
%set(gcf,'Position', [898 1011 697 422]);
set(gcf,'Position', [898 984 697 449]);

% set boundary conditions
ep(1)=0.01;
ep(2)=0.1;
u1=1;
u2=0;
x0=-1;
x1=4;
v0=0.5*(u2+u1);
w0=0.5*(u2-u1);

% parameters for plot
nx = 2000;
nt=3;
t(1)=0;
t(2)=2;
t(3)=4;
t(4)=4;

x=linspace(x0,x1,nx); 

subaxis(3,1,1,1,'MT',0.005,'MB',0.07,'MR',0.01,'ML',0.06,'P',0.03);
%subplot(3,1,1)
for i=1:nx
	u(i,1)=g(x(i),u1,u2);
end;
%plot(x,u(:,1))
plot([-1 -0.00001 0.0001, 4],[1 1 0 0],'Linewidth',1.1)
grid on
box on
axis([x0 x1 -0.18 1.18])
set(gca,'xtick',[-1 0 1 2 3 4]);
ylabel('u-axis','FontSize',12,'FontWeight','bold')
set(gca,'FontSize',12); 
%say=['\epsilon = ',num2str(ep)];
text(3.4,0.8,'t = 0','FontSize',12,'FontWeight','bold')

for it=2:nt

subaxis(3,1,1,it);
%subplot(3,1,it)
hold on
for iep=1:2
	for ix=1:nx
		q1 = 0.5*(x(ix)-u1*t(it))/sqrt(ep(iep)*t(it));
		q2 = 0.5*(x(ix)-u2*t(it))/sqrt(ep(iep)*t(it));
		q3 = (x(ix)-v0*t(it))*w0/ep(iep);
		k = erfc(q1)*exp(q3)/erfc(-q2);
		u(ix,it)=(u2+k*u1)/(1+k);
	end;
	if iep==1
		plot(x,u(:,it),'Linewidth',1.1);
	else
		plot(x,u(:,it),'--','Linewidth',1.1);
	end;
	end
	grid on
	box on
	if it==3
		xlabel('x-axis','FontSize',12,'FontWeight','bold')
		text(3.4,0.8,'t = 4','FontSize',12,'FontWeight','bold')
	else
		text(3.4,0.8,'t = 2','FontSize',12,'FontWeight','bold')
	end
	ylabel('u-axis','FontSize',12,'FontWeight','bold')
	legend(' \epsilon = 0.01',' \epsilon = 0.1',3)
	set(gca,'xtick',[-1 0 1 2 3 4]);
	axis([x0 x1 -0.18 1.18])

	% Set the fontsize to 12 for the plot 
	set(gca,'FontSize',12); 
	set(findobj(gcf,'tag','legend'),'FontSize',12,'FontWeight','bold');
end
hold off


function q=g(x,u1,u2)
if x < 0
	q=u1;
else
	q=u2;
end;