Octave: numerical solution of ODE

 

#  phase picture for the equation y'=2yx-2, 
#	where y=y(x) is the unknown function


## name some colors ...
blue="b";
red="r";

### more graphs in one picture ...
hold("on");

## the area of the plot ...
xx=4;
yy=3;
xlim([-xx/2,xx]);
ylim([-yy/2,yy]); 


## no axes around to be seen ...
axis("off");

### coordinate axes: dotted in black ...
dx=0.05;
t=-xx:dx:xx;
plot(t,0*t,".k");
plot(0*t,t,".k");

### the set {y'=0}, i.e. y=1/x
# (red color)
x=linspace(0.1,xx,30);
plot(x,1./x,red);
plot(-x,-1./x,red);

### the set {y"=0}, i.e. y = x/(1+2x^2)
# (blue color)
dd = @(x) x/(1+2*x*x);
x=linspace(-xx,xx,60);
plot(x,arrayfun(dd,x),blue);

### compute (numerically) the orbits 
h=0.01;
## p.s. ODR
function ydot=f(y,x)
  ydot=2*y*x-2;
endfunction
## initial condition are given in this cycle:
for ic = [1.75,2,0,-1,-2.2]
	## forward trajectory:
	x=0:h:xx;
	y=lsode("f",ic,x);
	plot(x,y,"k");
	## backward trajectory:
	x=0:-h:-xx;
	y=lsode("f",ic,x);
	plot(x,y,"k");
## end of the cycle ...
endfor

## save the plot as jpg ...
print -djpg obr2.jpg