
%
%   Script en Matlab para trazar la curvas solucion
%   de la funcion implicita
%
%   f(x,y)=x^3-y
%
%   usando el metodo de Newton.
% 
%   Rutinas a usar:
%
%    fun1.m       m-file donde se define la funcion en x
%    funder1.m    m-file donde se define la derivada de la funcion en x
%    newton.m    drive m-file que resuelve el cero de una funcion por medio
%                del Metodo de Newton
%
%    Metodologia:
%                Construimos una malla sobre la casa para visualizar las curvas.
%                Barremos para y de abajo hacia arriba, barremos en x de izquierda a derecha
%                considerando como punto inicial el punto de la malla.
%                Si Newton es exitoso, el punto resultado pertenece a la curva.
%
global y0

close all   % cerramos todas la ventanas

%
% Definimos la caja donde deseamos visualizar la curva
%

a=-5.;   % extremo izquierdo en x
b=5;     % extremo derecho en x
c=-5;    % extremo izquierdo en y
d=5;     % extremo derecho en y

%
% Definimos la dimension de la rejilla
%
gridx=25;   % numero de lineas verticales
gridy=55;   % numero de lineas horizontales
%
% Construimos los puntos igualmente espaciados en x y en y para  
% considerar la regilla sobre la caja de visualizacion
%
ab=linspace(a,b,gridx); 
cd=linspace(c,d,gridy);
%
% Asignamos el valor de los parametros
% para el Metodo de Newton
%
tolx=10^-4;    % Relative error
% epsilon=10^-8;  % Value to accept y0 as zero
nmax=25;   % Maximun number of iterations
%
% Inicializamos la matriz donde hemos de guardar los puntos  
% de la curva
%
pg=[];
%
% Barremos los puntos de la malla para identificar los puntos de
% la curva
%
for j=1:gridy,
   y0=cd(j);  
   for i=1:gridx,
      x0=ab(i);     % punto inicial para Newton
      [xs,niter,iflag] = newton('fun1','funder1',x0,nmax,tolx);
      if (iflag==5),
      %
      % Si converge el metodo, el punto [xs,y0] pertence a la curva
      %
         pg=[pg; xs y0]; % guardamos la informacion en una variable
      end
   end
end
%
% Graficamos los puntos que pertencen a la curva
%
plot(pg(:,1),pg(:,2),'.'); axis([a b c d]); grid;
xlabel('eje x'); ylabel('eje y'); title('curvas solucion');
