function [t,X]=odenullcline(TIME,xinit,A,I);
%%     Carol Lucas
%%     Subroutine called by AppendA_Dem
%%     Needs clode2 as calling program for ODE program

a11=A(1,1);a12=A(1,2);a21=A(2,1);a22=A(2,2);
T=trace(A)
D=det(A)
lambdas=eig(A)

exact=-inv(A)*I

[t,X] =ode45('clode2',TIME,xinit,[],A,I);

%% Much effort just to get the plotting ranges in the ball park

x1min=min(X(:,1))-10;x1min=min(x1min,exact(1)-10);
x1max=max(X(:,1))+10;x1max=max(x1max,exact(1)+10);
x2min=min(X(:,2))-10;x2min=min(x2min,exact(2)-10);
x2max=max(X(:,2))+10;x2max=max(x2max,exact(2)+10);
x1mintest=(-a12*x2min-I(1))/a11;
x1maxtest=(-a12*x2max-I(1))/a11;
x2mintest=(-a21*x1min-I(2))/a22;
x2maxtest=(-a21*x1max-I(2))/a22;
x1min=min(x1min,min(x1mintest,x1maxtest));
x1max=max(x1max,max(x1maxtest,x1mintest));
x2min=min(x2min,min(x2mintest,x2maxtest));
x2max=max(x2max,max(x2maxtest,x1mintest));
fmax=max(x1max,x2max);

x1ray=[x1min:x1max];x2ray=[x2min:x2max];

x1=(-a12*x2ray-I(1))/a11;
x2=(-a21*x1ray-I(2))/a22;

subplot(121);plot(t,X);xlabel('time');ylabel('X');
subplot(122);plot(X(:,1),X(:,2),x1,x2ray,x1ray,x2);hold on;
plot(exact(1),exact(2),'r*',X(1,1),X(1,2),'md');
xlabel('x1');ylabel('x2');title('nullclines');axis([x1min x1max x2min x2max]);
text(x1min,.9*x2max+.1*x2min,'  A = ');
text(.8*x1min+.2*x1max,.9*x2max+.1*x2min,num2str(A,0));
text(x1min,.1*x2max+.9*x2min,'  lambdas = ')
text(.6*x1min+.4*x1max,.1*x2max+.9*x2min,num2str(lambdas,2));