deff('yprim=f(t,y)', 'yprim=y^2-t')
t=0:0.1:5;
sol=ode(0,0,t,f);
plot2d(t, sol)
t=0:0.5:5; nt=length(t);
y=-3:3; ny=length(y);
ft=ones(nt, ny);
fy=feval(t,y,f);
xbasc(); champ(t,y,ft,fy)
plot2d(0:0.1:5, sol,frameflag=0)


function yprim= volterra(t, y)
yprim(1) = a*y(1) - b*y(1)*y(2);
yprim(2) = b*y(1)*y(2) - d*y(2);
endfunction  
a=3; b=1; d=2;

sol=ode([2;2], 0, 0:0.05:3, volterra);  
subplot(211); plot2d(0:0.05:3, sol(1,:)) 
subplot(212); plot2d(0:0.05:3, sol(2,:))
xbasc(); param3d(0:0.05:3, sol(1,:), sol(2,:))

xset('window',1); fchamp(volterra, 0, 0:0.2:4, 1:0.2:5)
plot2d(sol(1,:), sol(2,:), frameflag=0)             



function [f]=vpol(t, u)
f(1) = u(2);
f(2) = mu*(1-u(1)^2)*u(2) - u(1);
endfunction
 
 
mu=0.4;  
sol=ode([-2.5; 2.5], 0, 0:0.1:30, vpol);
xset('window', 0);
subplot(211); plot2d(0:0.1:30, sol(1,:))
subplot(212); plot2d(0:0.1:30, sol(2,:))
xset('window', 1);
plot(sol(1,:), sol(2,:))


deff('J=jac(t,y)', 'J=[a-b*y(2) -b*y(1); b*y(2) b*y(1)-d]')
t=0:0.1:5;
sol=ode([2; 2], 0, t, volterra, jac);
sol1=ode(2[; 2], 0, t, volterra);   
max(abs(sol1-sol))