It is an excelent job, however I found a bug. The function finds wrong intersection, if one of the curves is almost vertical. Please try this:
x1=[0.1 0.1+1e-17]
y1=[-1 1]
x2=[-1 1]
y2=[-1 1]
[x,y]=curveintersect(x1,y1,x2,y2)
plot(x1,y1,'k',x2,y2,'b:',x,y,'ro')
