% This is a simple ode integrator  
% it is set up to be easily modified to use for different dynamics
function  desolver

% plot ``nullclines"
figure(1)
plot([-2,2],[-2,2],'--',[-2,2],[2,-2],'--','linewidth',2)
hold on
 
tstep = .1;
t_start = -2;
t_end = 2;

for j = 1:21
 
 X0 = -2+0.2*(j-1);


%specify the output points
tspan = [t_start:tstep:t_end];

s0 = [X0];  %  format s0 = [U;V;...] as a column vector

[T,S] = ode23s(@deRHS,tspan, s0, odeset('maxstep',1));  
X = S(:,1); 
 
figure(1)
plot(T, X,-T,-X,'linewidth',2)

end

for j = 1:18
    X0 = -2;
    t_start = -2+0.2*(j-1);
    %specify the output points
tspan = [t_start:tstep:t_end];

s0 = [X0];  %  format s0 = [U;V;...] as a column vector

    [T,S] = ode23s(@deRHS,tspan, s0, odeset('maxstep',1));  
X = S(:,1); 
 
figure(1)
plot(T, X,-T,-X,'linewidth',2)

end

    
xlabel('t','fontsize',20)
ylabel('x','fontsize',20)
axis([-2 2 -2 2])
hold off


%the right hand side for ode simulation:
function s_prime=deRHS(t,s)
 
x = s(1);
 
 
s_prime = [x^2-t^2]';  
