clear

xbasc()

function res=v(x)
  res=0.5*x*x;
endfunction

function res1=syst1(x,y)
  e=y(3);
  res(3)=0.0;
  res(1)=y(2);
  res(2)=2*(v(x)-e)*y(1);
  res1=res';
endfunction

function res=psidr(e,xfin)
  xmin=-19.00;
  dx=0.001;
  y0=[0 1e-6 e]';
  x=xmin:dx:xfin;
  y=ode(y0,xmin,x,syst1);
//plot(x,y);
  res=y(2,size(y));
endfunction


function res=trace(e)
  xmin=-5.00;
  xmax=5.00;
  dx=0.01;
  y0=[0 1e-6 e]';
  x=xmin:dx:xmax;
  y=ode(y0,xmin,x,syst1);
  plot(x,y);
  res=y(2,size(y));
endfunction


function res=psiga(e,xfin)
  xmax=19.00;
  y0=[0 1e-6 e]';
  dx=0.001;
  x=xmax:-dx:xfin;
  y=ode(y0,xmax,x,syst1);
//plot(x,y);
  res=y(2,size(y));
endfunction

function res=fz(x)
  e=x;
//  disp(e)
  xd=4.23;
  dx=0.001;
  res=(psidr(e,xd+dx)-psidr(e,xd-dx))/psiga(e,xd)-(psiga(e,xd+dx)-psiga(e,xd-dx))/psidr(e,xd);
//  res=res/psidr(e,xd);
//  disp(res);
endfunction

i=0

for e=0:0.1:5
  i=i+1;
  f(i,2)=log(abs(fz(e)));
  f(i,1)=e;
  plot(f(:,1),f(:,2))
end

//eps=1;
//a=0.5;
//b=1.7;
//
//while(eps>1e-6)
//  c=(a+b)/2
//  eps=b-a;
//  if(fz(c)*fz(a)<0) 
//    b=c;
//   else
//    a=c;
//  end
//end

for e=0:0.3:5
  es=fsolve(e,fz)
  xbasc()
  trace(es)
  pause
end

