% plotxytf.m 

function p=plotxytf(x1func,x2func,t0,tf,linetype,plotxlabel) 
% Used by SYS3PLOT.M with PPLANE2F.M to plot user-derives solution(s) to an 
% ODE system overlayed on a phase plane plot. 

plotN=50; 
x1min=0; x2min=x1min; 
x1max=1; x2max=x1max; 
t=t0:(tf-t0)/(plotN-1):tf; 
ttemp=t; 
t=tf; 
if eval(x1func)>-inf 
	x1min=eval(x1func); 
	if eval(x1func)<inf 
		x1max=eval(x1func); 
		x1funcplot=[]; 
		t=[]; 
		for j=1:length(ttemp) 
			t=ttemp(j); 
			x1funcplot=[x1funcplot,eval(x1func)]; 
			if min(x1funcplot)>-inf 
				if x1min>min(x1funcplot) 
					x1min=min(x1funcplot); 
				end 
			end 
			if max(x1funcplot)<inf 
				if x1max<max(x1funcplot) 
					x1max=max(x1funcplot); 
				end 
			end 
		end 
	else 
		disp('Your function for x1 has a + inf result somewhere') 
		break 
	end 
else 
	disp('Your function for x1 has a - infinity result somewhere') 
	break 
end 

if eval(x2func)>-inf 
	x2min=eval(x2func); 
	if eval(x2func)<inf 
		x2max=eval(x2func); 
		x2funcplot=[]; 
		t=[]; 
		for j=1:length(ttemp) 
			t=ttemp(j); 
			x2funcplot=[x2funcplot,eval(x2func)]; 
			if min(x2funcplot)>-inf 
				if x2min>min(x2funcplot) 
					x2min=min(x2funcplot); 
				end 
			end 
			if max(x2funcplot)<inf 
				if x2max<max(x2funcplot) 
					x2max=max(x2funcplot); 
				end 
			end 
		end 
	else 
		disp('Your function for x2 has a + inf result somewhere') 
		break 
	end 
else 
	disp('Your function for x2 has a - infinity result somewhere') 
	break 
end 
ax=[x1min,x1max,x2min,x2max]; 
plot(x1funcplot,x2funcplot,linetype) 
axis(ax) 
xlabel(plotxlabel)