function cilindro
% Progama basico
%
% Este programa calcula el estado básico de un problema de Benard-Marangoni 
% en un dominio anular. 
%
% --- Punto 0 ----. Fijamos los parametros
n = 13;
m = 11;
dif = 2;  	 % Diferencia de temperatura entre las paredes
biot = 1.25; % Parametro de Biot 
tinicio = 3; 
tinicio2 = 2; 
diri = 0;				% aislante 
long = 0.002;
C = datos(n,m);

% --- Construimos el bucle de aproximación
[ur0,uz0,t0,p0] = estbas(n,m,long,biot,tinicio2,dif,0,C);
%
% Aproximo hasta tinicio
%
for deltat=linspace(tinicio2,tinicio,5)
	fprintf('-------------------------------\n')
	fprintf('calculando para deltat=%5.3f  \n',deltat)
	fprintf('-------------------------------\n')
	[ur0,uz0,t0,p0] = estbas(n,m,long,biot,deltat,dif,2,C);
end

save ('datos.fluidos.mat','-v6')

function [ur0,uz0,t0,p0] = estbas(n,m,long,biot,deltat,dif,preg,C)
% Programa estbas
% Este programa calcula el estado basico de un problema de colocación en un cilindro con
% condiciones de contorno Dirichlet o Neumann en las paredes. 
% La condición de paro esta sobre el residuo.
% El programa chequea los datos que solo dependen de la dimensión por si ya están grabados 
% para no calcularlos.
% --- Entradas --- 
% n		-> número de polinomios en la dirección radial   -> n+1 puntos de colocación en r
% m		-> número de polinomios en la dirección vertical -> m+1 puntos de colocación en z
% long 	-> áltura del cilindro.
% biot	-> Número de biot
% deltat	-> diferencia de temperatura desde la pared caliente a la fría
% dif 	-> Incremento de temperatura en la pared fría
% preg	-> Detalle técnico del algoritmo 0 para iterar desde aprox2, 1 desde el conductivo, 
%				otra cosa desde la anterior.
% diri	-> Si diri = 1, toma condiciones dirichlet. Si es distinto toma condiciones Neumann.
% 				Por defecto es 1. 
% [ur0,uz0,t0,p0] = definitivo.
%
% --- Rutinas ---
% Derevol para calcular las derivadas de las funciones
% Derpoly para calcular los polinomios de chebychev
% aprox2  para calcular una aproximación de la iteración.
% --- Rutinas asociadas ---
% delcrit 	-> Halla el deltat critico. Toma dos salidas
% --- Programa ---
%
maximo = max(n,m);
hi = 1;
diri = 0;
prfin = 0;
r = C{1};
z = C{2} ;
cheby = C{3};
dcheby = C{4};
ddcheby = C{5};
gcheby = C{6};
dgcheby = C{7};
ddgcheby = C{8};
%
% Ahora vamos a definir los parametros particulares de cada problema
%
aa = 0.01;    		  		    		    	% radio interior              
d = 0.02;										% radio de la corona	
a = 2*long/d;									% Aspect-Ratio
alpha = .001;									% thermal expansion
gamrho = 8.E-7/885.	;						% surface tension/density
nu = 3.22E-6	;								% viscosity
kappa = 0.7E-7	;								% thermal cond. diff
g = 9.81	;										% Constante de la gravedad
kn = kappa*nu;									% Cálculo técnico
T0 = 10;											% Temperatura ambiental
Tmax = T0+deltat;								% Temperatura en la pared interior 
Tmin = T0+deltat-dif;			   		% Temperatura en la pared exterior
sigma = nu/kappa;   							% Número de Prandtl
isigma = 1/sigma;		           			% Inversa del número de Prandtl
Ra = g.*deltat.*alpha.*long^3/kn; 		% Número de Rayleigh	
Mar = gamrho*deltat.*long/kn;				% Constante de Marangoni	
Bond = Ra/Mar;									% Número de Bond
b = long^3.*g/kn;								% Constante de cálculo
% Parámetros técnicos del problema
N = 2.*long./(2.*aa+d*(1+r));
qr = N(ones(1,maximo+1),:);
a2 = aa;
ja = 0; 											% Dato técnico
% Ahora definimos desde donde empezamos a iterar. 
if preg==0 		% aproximación lineal
   save (['datos_para_aprox_con_n=',num2str(n),'_m=',num2str(m)])
   [ur0,uz0,t0,p0] = aprox2(n,m);
else						% Estado anterior
   load datos1
end
% Predimensionalización de la matriz de datos y de los vectores del sistema.
ac = zeros(4*(n+1)*(m+1));
bc = zeros(4*(n+1)*(m+1),1);
%
%%%%%%%%%%%%%%%%%%%%%%% Matrices de la perturbación %%%%%%%%%%%%%%%%%%%%%%%%
%
% Construcción de la parte comun de la matriz que no cambia en las iteraciones. Sobre todo 
% se va a calcular la parte de las condiciones de contorno.
% Primero las ecuaciones de NS
ja = ja+(n-1)*(m-1);
ja = ja+(n-1)*(m-1);
% Tercera ecuación. Tambien cambia
ja = ja+(n-1)*(m-1);  
% Cuarta ecuación. Divergencia, ya no cambia
for i=2:m
   for J=1:m+1
      sumur1 = qr(1:n+1,2:n)'.*cheby(J,i).*gcheby(:,2:n)';
      sumur2 = a*cheby(J,i).*dgcheby(:,2:n)';
      ac(ja+(n-1)*(i-2)+(1:n-1),(J-1)*(n+1)+(1:n+1)) = sumur1+sumur2;
      sumuz = 2*dcheby(J,i).*gcheby(:,2:n)';    
      ac(ja+(n-1)*(i-2)+(1:n-1),(n+1)*(m+1)+(J-1)*(n+1)+(1:n+1)) = sumuz;
   end
end
ja=ja+(n-1)*(m-1);   
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Condiciones de contorno %
%%%%%%%%%%%%%%%%%%%%%%%%%%%
% En este caso se consideran casi todos
% Condición en r=-1 todos los zetas
% c.c. ur=0
for i=1:m+1
   for J=1:m+1
      ac(ja+i,(J-1)*(n+1)+(1:n+1))=cheby(J,i).*gcheby(1:n+1,1)';
   end
end
ja=ja+m+1;
% c.c uz=0
for i=1:m+1
   for J=1:m+1
      ac(ja+i,(n+1)*(m+1)+(J-1)*(n+1)+(1:n+1))=cheby(J,i).*gcheby(1:n+1,1)';
   end
end
ja=ja+m+1;
% Bucle de elección segun diri, tomando o dirichlet o neumann. por \partial_r T = 0
% c.c. T=1
for i=1:m+1
   for J=1:m+1
      if diri==1
         ac(ja+i,2*(n+1)*(m+1)+(J-1)*(n+1)+(1:n+1))=cheby(J,i).*gcheby(1:n+1,1)';
      else
         ac(ja+i,2*(n+1)*(m+1)+(J-1)*(n+1)+(1:n+1))=cheby(J,i).*dgcheby(1:n+1,1)';
      end
   end
end
ja=ja+m+1;
%   c.c. de Orszag sobre la presión.
ja=ja+m;
% Condiciones sobre r=1 (todos los z)
% c.c. ur=0
for i=1:m+1
   for J=1:m+1
      ac(ja+i,(J-1)*(n+1)+(1:n+1))=cheby(J,i).*gcheby(1:n+1,n+1)';
   end
end
ja=ja+m+1;   
% c.c uz=0
for i=1:m+1
   for J=1:m+1
      ac(ja+i,(n+1)*(m+1)+(J-1)*(n+1)+(1:n+1))=cheby(J,i).*gcheby(1:n+1,n+1)';
   end
end
ja=ja+m+1;
% Elijo segun marque diri para la temperatura
for i=1:m+1
   for J=1:m+1
      if diri==1
         ac(ja+i,2*(n+1)*(m+1)+(J-1)*(n+1)+(1:n+1)) = cheby(J,i).*gcheby(1:n+1,n+1)';
      else
         ac(ja+i,2*(n+1)*(m+1)+(J-1)*(n+1)+(1:n+1)) = cheby(J,i).*dgcheby(1:n+1,n+1)';
      end
   end
end
ja = ja+m+1;
%   c.c. de Orszag sobre la presión.
ja = ja+m;  
% Condiciones para la capa z=-1 (todos los r menos los de los bordes)  
% ur=0
for J=1:m+1 
   ac(ja+(1:n-1),(J-1)*(n+1)+(1:n+1))=cheby(J,1).*gcheby(:,2:n)';            
end
ja=ja+n-1;
%uz=0
for J=1:m+1
   ac(ja+(1:n-1),(n+1)*(m+1)+(J-1)*(n+1)+(1:n+1))=cheby(J,1).*gcheby(:,2:n)';
end
ja=ja+n-1;
% Decrecimiento lineal en la frontera
for J=1:m+1
   ac(ja+(1:n-1),2*(n+1)*(m+1)+(J-1)*(n+1)+(1:n+1))=cheby(J,1).*gcheby(:,2:n)';
end
ja=ja+n-1;
% Condición divergencia nula
for J=1:m+1
   sumur1=qr(:,2:n)'.*cheby(J,1).*gcheby(:,2:n)';
   sumur2=a*cheby(J,1).*dgcheby(:,2:n)';
   ac(ja+(1:n-1),(J-1)*(n+1)+(1:n+1))= sumur2+sumur1;
   sumuz=2*dcheby(J,1).*gcheby(:,2:n)';    
   ac(ja+(1:n-1),(n+1)*(m+1)+(J-1)*(n+1)+(1:n+1))=sumuz;  
end
ja=ja+n-1;
% Condiciones para la capa z=1 (todos los r menos los de los bordes)
% 2*\partial_z ur+M.*a.*partial_r T=0
for J=1:m+1
   ac(ja+(1:n-1),(J-1)*(n+1)+(1:n+1))=2*dcheby(J,m+1).*gcheby(:,2:n)';
   ac(ja+(1:n-1),2*(n+1)*(m+1)+(J-1)*(n+1)+(1:n+1))=Mar.*a.*cheby(J,m+1).*dgcheby(:,2:n)';
end
ja=ja+n-1;   
% uz=0
for J=1:m+1
   ac(ja+(1:n-1),(n+1)*(m+1)+(J-1)*(n+1)+(1:n+1)) = cheby(J,m+1).*gcheby(:,2:n)';
end
ja=ja+n-1;
% Condición de biot   
for J=1:m+1
   sumT1 = 2*dcheby(J,m+1).*gcheby(:,2:n)';
   sumT2 = biot.*cheby(J,m+1).*gcheby(:,2:n)';
   ac(ja+(1:n-1),2*(n+1)*(m+1)+(J-1)*(n+1)+(1:n+1)) = sumT1+sumT2;
end
% Fin de la construcción de la parte comun de la matriz
% Ahora construimos el bucle principal
% Inicio las normas
n1 = 1;
n2 = 0;
n3 = 0;
n4 = 0;
paro = 20;
while (hi<paro)& (n1+n2+n3+n4>1e-7)
   % Empezamos construyendo las evaluaciones de las derivadas.
   [drur0,dzur0,drrur0,dzzur0,Trur0,Tzur0] = derevol(ur0,C);
   [druz0,dzuz0,drruz0,dzzuz0,Truz0,Tzuz0] = derevol(uz0,C);
   [drt0,dzt0,drrt0,dzzt0,Trt0,Tzt0] = derevol(t0,C);
   [drp0,dzp0] = derevol(p0,C);
   Tur0 = zeros(n+1,n+1,m+1);
   Tuz0 = zeros(n+1,n+1,m+1);
   for i=1:m+1
      auxr=ur0(i,:);auxz=uz0(i,:);
      Tur0(:,:,i) = auxr(ones(maximo+1,1),:);
      Tuz0(:,:,i) = auxz(ones(maximo+1,1),:);
   end
   % Primera ecuación.
   ja = 0;
   for i=2:m
      for J=1:m+1
         sump = -a*cheby(J,i).*dgcheby(:,2:n)';
         ac(ja+(n-1)*(i-2)+(1:n-1),3*(n+1)*(m+1)+(J-1)*(n+1)+(1:n+1)) = sump;
         sumur1 = a*a*cheby(J,i).*ddgcheby(:,2:n)';
         sumur2 = a*qr(1:n+1,2:n)'*cheby(J,i).*dgcheby(:,2:n)';
         sumur3 = 4*ddcheby(J,i).*gcheby(:,2:n)';
         sumur4 = -qr(1:n+1,2:n)'.^2.*cheby(J,i).*gcheby(:,2:n)';
         if prfin==1							% Si consideramos Pr finito
            nolin1 = -isigma*a*Tur0(1:n+1,2:n,i)'.*cheby(J,i).*dgcheby(:,2:n)';
            nolin2 = -isigma*a*Trur0(1:n+1,2:n,i)'.*cheby(J,i).*gcheby(:,2:n)';
            nolin3 = -2*isigma*Tuz0(1:n+1,2:n,i)'.*dcheby(J,i).*gcheby(:,2:n)';
            nolin4 = -2*isigma*Tzur0(:,2:n,i)'.*cheby(J,i).*gcheby(:,2:n)';        
            ac(ja+(n-1)*(i-2)+(1:n-1),(J-1)*(n+1)+(1:n+1)) = sumur1+sumur2+sumur3+...
               sumur4+nolin1+nolin2+nolin3;
            ac(ja+(n-1)*(i-2)+(1:n-1),(n+1)*(m+1)+(J-1)*(n+1)+(1:n+1)) = nolin4;
         else
            ac(ja+(n-1)*(i-2)+(1:n-1),(J-1)*(n+1)+(1:n+1)) = sumur1+sumur2+sumur3+sumur4;
         end
      end
   end
   ja = ja+(n-1)*(m-1);
   % Creación del termino independiente.
   ind = (a*drp0(2:m,2:n)-a.*a*drrur0(2:m,2:n)-qr(2:m,2:n).*a.*drur0(2:m,2:n)...
      -4*dzzur0(2:m,2:n)+qr(2:m,2:n).*qr(2:m,2:n).*ur0(2:m,2:n))';

   if prfin==1
      ind = ind+(isigma*ur0(2:m,2:n).*a.*drur0(2:m,2:n)+...
         2*isigma*uz0(2:m,2:n).*druz0(2:m,2:n))';
   end
   bc(1:ja) = ind(:);
   % Segunda ecuacion. 
   for i=2:m
      for J=1:m+1
         sump = -2*dcheby(J,i).*gcheby(:,2:n)';
         ac(ja+(n-1)*(i-2)+(1:n-1),3*(n+1)*(m+1)+(J-1)*(n+1)+(1:n+1)) = sump;
         sumuz1 = a*a*cheby(J,i).*ddgcheby(:,2:n)';
         sumuz2 = a*qr(1:n+1,2:n)'.*cheby(J,i).*dgcheby(:,2:n)';
         sumuz3 = 4*ddcheby(J,i).*gcheby(:,2:n)'; 
         sumt = Ra.*cheby(J,i).*gcheby(:,2:n)';           
         ac(ja+(n-1)*(i-2)+(1:n-1),2*(n+1)*(m+1)+(J-1)*(n+1)+(1:n+1)) = sumt; 
         if prfin==1
            nolin1=-isigma*a*Tur0(1:n+1,2:n,i)'.*cheby(J,i).*dgcheby(:,2:n)';
            nolin2=-isigma*a*Truz0(1:n+1,2:n,i)'.*cheby(J,i).*gcheby(:,2:n)';
            nolin3=-2*isigma*Tuz0(1:n+1,2:n,i)'.*dcheby(J,i).*gcheby(:,2:n)';
            nolin4=-2*isigma*Tzuz0(:,2:n,i)'.*cheby(J,i).*gcheby(:,2:n)';
            ac(ja+(n-1)*(i-2)+(1:n-1),(J-1)*(n+1)+(1:n+1))=nolin2;
            ac(ja+(n-1)*(i-2)+(1:n-1),(n+1)*(m+1)+(J-1)*(n+1)+(1:n+1))=sumuz1+sumuz2+sumuz3+...
               nolin1+nolin3+nolin4;
         else
            ac(ja+(n-1)*(i-2)+(1:n-1),(n+1)*(m+1)+(J-1)*(n+1)+(1:n+1))=sumuz1+sumuz2+sumuz3;
         end      
      end
   end
   ind = (2*dzp0(2:m,2:n)-a*a*drruz0(2:m,2:n)-a*qr(2:m,2:n).*druz0(2:m,2:n)-...
      4.*dzzuz0(2:m,2:n)-Ra.*t0(2:m,2:n))';
   if prfin==1
      ind = ind+(isigma*a*ur0(2:m,2:n).*druz0(2:m,2:n)+2*isigma*uz0(2:m,2:n).*dzuz0(2:m,2:n))';
   end
   bc(ja+(1:(n-1)*(m-1)))=b+ind(:);
   ja=ja+(n-1)*(m-1);
   %tercera ecuación:calor
   for i=2:m
      for J=1:m+1
         nolin1 = -a*Trt0(:,2:n,i)'.*cheby(J,i).*gcheby(:,2:n)';
         nolin2 = -a*Tur0(1:n+1,2:n,i)'.*cheby(J,i).*dgcheby(:,2:n)';
         nolin3 = -2*Tuz0(1:n+1,2:n,i)'.*dcheby(J,i).*gcheby(:,2:n)';
         nolin4 = -2*Tzt0(:,2:n,i)'.*cheby(J,i).*gcheby(:,2:n)';
         ac(ja+(n-1)*(i-2)+(1:n-1),(J-1)*(n+1)+(1:n+1)) = nolin1;
         ac(ja+(n-1)*(i-2)+(1:n-1),(n+1)*(m+1)+(J-1)*(n+1)+(1:n+1)) = nolin4;
         sumt1 = a*a.*cheby(J,i).*ddgcheby(:,2:n)';
         sumt2 = a.*qr(1:n+1,2:n)'.*cheby(J,i).*dgcheby(:,2:n)';
         sumt3 = 4.*ddcheby(J,i).*gcheby(:,2:n)';
         ac(ja+(n-1)*(i-2)+(1:n-1),2*(n+1)*(m+1)+(J-1)*(n+1)+(1:n+1)) = nolin2+nolin3...
            +sumt1+sumt2+sumt3;
      end
   end
   ind = (a*drt0(2:m,2:n).*ur0(2:m,2:n)+2*uz0(2:m,2:n).*dzt0(2:m,2:n)-a.*a.*drrt0(2:m,2:n)+...
      -a*qr(2:m,2:n).*drt0(2:m,2:n)-4.*dzzt0(2:m,2:n))';
   bc(ja+(1:(n-1)*(m-1))) = ind(:);
   ja = ja+(n-1)*(m-1);  
   % cuarta ecuación: divergencia. No la calculamos porque no cambia
   ind = (-a.*drur0(2:m,2:n)-qr(2:m,2:n).*ur0(2:m,2:n)-2*dzuz0(2:m,2:n))';
   bc(ja+(1:(n-1)*(m-1)))=ind(:);
   ja=ja+(n-1)*(m-1);   
   %%%%%%%%%%%%%%%%%%%%%%%%%%%
   % Condiciones de contorno %%%%
   %%%%%%%%%%%%%%%%%%%%%%%%%%%
   % Condición en r=-1 todos los zetas
   % c.c. ur=0
   bc(ja+(1:m+1))=-ur0(:,1);
   ja=ja+m+1;
   % c.c uz=0
   bc(ja+(1:m+1))=-uz0(:,1);
   ja=ja+m+1;
   % Bucle de elección segun diri, tomando o dirichlet o neumann. por \partial_r T = 0
   % c.c. T=1
   if diri==1
      bc(ja+(1:m+1))=1-t0(:,1);
   else
      bc(ja+(1:m+1))=-drt0(:,1);
   end
   ja=ja+m+1;
   %   c.c. de Orszag sobre la presión.
   for i=1:m
      for J=1:m+1
         sump = -a*cheby(J,i).*dgcheby(:,1)';
         ac(ja+i,3*(n+1)*(m+1)+(J-1)*(n+1)+(1:n+1))=sump;
         sumur1 = a*a*cheby(J,i).*ddgcheby(:,1)';
         sumur2 = a*qr(1:n+1,1)'*cheby(J,i).*dgcheby(:,1)';
         sumur3 = 4*ddcheby(J,i).*gcheby(:,1)';
         sumur4 = -qr(:,1)'.*qr(:,1)'.*cheby(J,i).*gcheby(:,1)';
         if prfin==1
            nolin1 = -isigma*a*Tur0(1:n+1,1,i)'.*cheby(J,i).*dgcheby(:,1)';
            nolin2 = -isigma*a*Trur0(1:n+1,1,i)'.*cheby(J,i).*gcheby(:,1)'; 
            nolin3 = -2*isigma*Tuz0(1:n+1,1,i)'.*dcheby(J,i).*gcheby(:,1)';
            nolin4 = -2*isigma*Tzur0(:,1,i)'.*cheby(J,i).*gcheby(:,1)';
            ac(ja+i,(J-1)*(n+1)+(1:n+1))=sumur1+sumur2+sumur3+sumur4+...
               nolin1+nolin2+nolin4+nolin3;
         else
            ac(ja+i,(J-1)*(n+1)+(1:n+1))=sumur1+sumur2+sumur3+sumur4;
         end   
      end
   end
   clear ind;
   % No se pone traspuesta porque no cambia, es un vector.
   ind = a*drp0(1:m,1)-a.*a*drrur0(1:m,1)-qr(1:m,1).*a.*drur0(1:m,1)-4*dzzur0(1:m,1)+...
      qr(1:m,1).*qr(1:m,1).*ur0(1:m,1);
   if prfin==1
      ind = ind+isigma*a*ur0(1:m,1).*drur0(1:m,1)+2*isigma*uz0(1:m,1).*dzur0(1:m,1);
   end
   bc(ja+(1:m))=ind(:);clear ind,
   ja=ja+m;
   % Condiciones sobre r=1 (todos los z)
   % c.c. ur=0
   bc(ja+(1:m+1))=-ur0(:,n+1);
   ja=ja+m+1;   
   % c.c uz=0
   bc(ja+(1:m+1))=-uz0(:,n+1);
   ja=ja+m+1;
   % Elijo segun marque diri para la temperatura
   if diri==1
      bc(ja+(1:m+1))=(Tmin-T0)./deltat-t0(:,n+1);
   else
      bc(ja+(1:m+1))=-drt0(:,n+1);
   end
   ja=ja+m+1;
   %   c.c. de Orszag sobre la presión.
   for i=1:m
      for J=1:m+1
         sump = -a*cheby(J,i).*dgcheby(:,n+1)';
         ac(ja+i,3*(n+1)*(m+1)+(J-1)*(n+1)+(1:n+1))=sump;
         sumur1 = a*a*cheby(J,i).*ddgcheby(:,n+1)';
         sumur2 = a*qr(1:n+1,n+1)'*cheby(J,i).*dgcheby(:,n+1)';
         sumur3 = 4*ddcheby(J,i).*gcheby(:,n+1)';
         sumur4 = -qr(:,n+1)'.*qr(:,n+1)'.*cheby(J,i).*gcheby(:,n+1)';
         if prfin==1
            nolin1 = -isigma*a*Tur0(1:n+1,n+1,i)'.*cheby(J,i).*dgcheby(:,n+1)';
            nolin2 = -isigma*a*Trur0(1:n+1,n+1,i)'.*cheby(J,i).*gcheby(:,n+1)'; 
            nolin3 = -2*isigma*Tuz0(1:n+1,n+1,i)'.*dcheby(J,i).*gcheby(:,n+1)';
            nolin4 = -2*isigma*Tzur0(:,n+1,i)'.*cheby(J,i).*gcheby(:,n+1)';
            ac(ja+i,(J-1)*(n+1)+(1:n+1))=sumur1+sumur2+sumur3+sumur4+...
               nolin1+nolin2+nolin4+nolin3;
         else
            ac(ja+i,(J-1)*(n+1)+(1:n+1))=sumur1+sumur2+sumur3+sumur4;
         end   
      end
   end
   clear ind;
   ind = a*drp0(1:m,n+1)-a.*a*drrur0(1:m,n+1)-qr(1:m,n+1).*a.*drur0(1:m,n+1)-...
      4*dzzur0(1:m,n+1)+qr(1:m,n+1).*qr(1:m,n+1).*ur0(1:m,n+1);
   if prfin==1
      ind = ind+isigma*a*ur0(1:m,n+1).*drur0(1:m,n+1)+2*isigma*uz0(1:m,n+1).*dzur0(1:m,n+1);
   end
   bc(ja+(1:m)) = ind(:);clear ind
   for J=1:m+1      
      ac(ja+4,(J-1)*(n+1)+(1:n+1))=0;
      ac(ja+4,(n+1)*(m+1)+(J-1)*(n+1)+(1:n+1))=0;
      ac(ja+4,3*(n+1)*(m+1)+(J-1)*(n+1)+(1:n+1))=cheby(J,4).*gcheby(:,n+1)';
   end   
   bc(ja+4) = -p0(4,1);					% Ojo, esto es arbitrario.
   ja = ja+m;  
   % Condiciones para la capa z=-1 (todos los r menos los de los bordes)  
   % ur=0
   bc(ja+(1:n-1))=-ur0(1,2:n);
   ja=ja+n-1;
   %uz=0
   bc(ja+(1:n-1))=-uz0(1,2:n);
   ja=ja+n-1;
   % Decrecimiento lineal en la frontera
   bc(ja+(1:n-1))=(Tmin-Tmax)*long/d/deltat*(1./N(2:n)-a2/long)+1-t0(1,2:n);
   ja=ja+n-1;
   % Condición divergencia nula
   bc(ja+(1:n-1))=-drur0(1,2:n)*a-2*dzuz0(1,2:n)-qr(1,2:n).*ur0(1,2:n);
   ja=ja+n-1;
   %
   % Condiciones para la capa z=1 (todos los r menos los de los bordes)
   %
   % 2*\partial_z ur+M.*a.*partial_r T=0
   bc(ja+(1:n-1))=-2*dzur0(m+1,2:n)'-Mar.*a.*drt0(m+1,2:n)';
   ja=ja+n-1;   
   % uz=0
   bc(ja+(1:n-1))=-uz0(m+1,2:n)';
   ja=ja+n-1;
   % Condición de biot   
   bc(ja+(1:n-1))=-2*dzt0(m+1,2:n)'-biot.*t0(m+1,2:n)';
   ja=ja+n-1;
   
   % Condición Orszag
   for J=1:m+1
      sump = -2*dcheby(J,m+1).*gcheby(:,1:n+1)';
      ac(ja+(1:n+1),3*(n+1)*(m+1)+(J-1)*(n+1)+(1:n+1))=sump;
      sumuz1 = a*a*cheby(J,m+1).*ddgcheby(:,1:n+1)'; 
      sumuz2 = a*qr(1:n+1,1:n+1)'.*cheby(J,m+1).*dgcheby(:,1:n+1)';
      sumuz3 = 4*ddcheby(J,m+1).*gcheby(:,1:n+1)';
      if prfin==1
         nolin1 = -isigma*a*Tur0(1:n+1,:,i)'.*cheby(J,m+1).*dgcheby(:,:)';
         nolin2 = -isigma*a*Truz0(1:n+1,:,i)'.*cheby(J,i).*gcheby(:,:)';
         nolin3 = -2*isigma*Tuz0(1:n+1,:,i)'.*dcheby(J,i).*gcheby(:,:)';
         nolin4 = -2*isigma*Tzuz0(:,:,i)'.*cheby(J,i).*gcheby(:,:)';
         ac(ja+(1:n+1),(J-1)*(n+1)+(1:n+1))=nolin2;
         ac(ja+(1:n+1),(n+1)*(m+1)+(J-1)*(n+1)+(1:n+1))=sumuz1+sumuz2+sumuz3...
            +nolin1+nolin3+nolin4;
      else
         ac(ja+(1:n+1),(n+1)*(m+1)+(J-1)*(n+1)+(1:n+1))=sumuz1+sumuz2+sumuz3;
      end   
      sumt=Ra.*cheby(J,m+1).*gcheby(:,1:n+1)';    
      ac(ja+(1:n+1),2*(n+1)*(m+1)+(J-1)*(n+1)+(1:n+1))=sumt;   
   end
   ind = 2*dzp0(m+1,1:n+1)'-a*a*drruz0(m+1,1:n+1)'-a*qr(m+1,1:n+1)'.*druz0(m+1,1:n+1)'...
      -4.*dzzuz0(m+1,1:n+1)'-Ra.*t0(m+1,1:n+1)';
   if prfin==1
      ind = ind+isigma*a*ur0(m+1,:)'.*druz0(m+1,:)'+2*isigma*uz0(m+1,:)'.*dzuz0(m+1,:)';
   end
   bc(ja+(1:n+1))=b+ind(:);
   % Fin de la construcción de las matrices. Resolvemos y calculamos las perturbaciones
   solucion=ac\bc;
   coefur = reshape(solucion(1:(n+1)*(m+1)),n+1,m+1)';
   coefuz = reshape(solucion((n+1)*(m+1)+1:2*(n+1)*(m+1)),n+1,m+1)';
   coeft = reshape(solucion(2*(n+1)*(m+1)+1:3*(n+1)*(m+1)),n+1,m+1)';
   coefp = reshape(solucion(3*(n+1)*(m+1)+1:4*(n+1)*(m+1)),n+1,m+1)';
   % Predimensionalización de las matrices de datos
   
   Ur = (cheby')*coefur*gcheby;
   Uz = (cheby')*coefuz*gcheby;
   T = (cheby')*coeft*gcheby;
   P = (cheby')*coefp*gcheby;
   %
   % Las normas calculan el tamaño de estas perturbaciones y solo de las perturbaciones
   %
   n1=norm(Ur,2);
   n2=norm(Uz,2);
   n3=norm(T,2);
   n4=norm(P,2);
	fprintf('paso %i: n1=%e, n2=%e n3=%e n4=%e \n',hi,n1,n2,n3,n4)
   %if (n1>1e1) %| (hi>10)
	%   fprintf('paso %i: n1=%e, n2=%e n3=%e n4=%e \n',hi,n1,n2,n3,n4)
   %end
   hi=hi+1;
   ur0=Ur+ur0;
   uz0=Uz+uz0;
   t0=T+t0;
   p0=P+p0;
end	
% Guardamos las soluciones que hemos calculado si converge. 
if hi<paro
   save datos1.mat ur0 uz0 t0 p0
end
%
% ---- Principio function aprox2
%
function [ur0,uz0,t0,p0] = aprox2(n,m)
% Programa definitivo
% Este programa calcula el autovalor máximo de un problema de colocación en un cilindro con
% condiciones de contorno Neumann. La condición de paro esta sobre el residuo
% --- Entradas --- 
% n		-> número de polinomios en la dirección radial   -> n+1 puntos de colocación en r
% m		-> número de polinomios en la dirección vertical -> m+1 puntos de colocación en z
% long 	-> áltura del cilindro.
% biot	-> Número de biot
% deltat	-> diferencia de temperatura desde la pared caliente a la fría
% dif 	-> Incremento de temperatura en la pared fría
% diri		-> Si diri = 1, toma condiciones dirichlet. Si es distinto toma condiciones Neumann.
% 				Por defecto es 1. 
% --- Salidas ---
% [ur0,uz0,t0,p0] 
% Donde estos son los valores de los campos negando las nolinealidades
% --- Rutinas ---
% Ninguna
% --- Programa ---
% Cargamos todos los datos técnicos. 
load (['datos_para_aprox_con_n=',num2str(n),'_m=',num2str(m)])
qr = N(ones(1,maximo+1),:);
% Predimensionalización de las matrices
ac = zeros(4*(n+1)*(m+1));
bc = zeros(size(ac(:,1)));
ja=0;
% Primera ecuación
for i=2:m
   for j=1:m+1
      sump = -a*cheby(j,i).*dgcheby(:,2:n)';
      ac(ja+(n-1)*(i-2)+(1:n-1),3*(n+1)*(m+1)+(j-1)*(n+1)+(1:n+1)) = sump;
      sumur1 = a*a*cheby(j,i).*ddgcheby(:,2:n)';
      sumur2 = a*qr(:,2:n)'.*cheby(j,i).*dgcheby(:,2:n)';
      sumur3 = 4*ddcheby(j,i).*gcheby(:,2:n)';
      sumur4 = -(qr(:,2:n)').^2.*cheby(j,i).*gcheby(:,2:n)';
      ac(ja+(n-1)*(i-2)+(1:n-1),(j-1)*(n+1)+(1:n+1)) = sumur1+sumur2+sumur3+sumur4;
   end
end
ja=ja+(n-1)*(m-1);
% Segunda ecuación.
for i=2:m
   for j=1:m+1
      sump = -2*dcheby(j,i).*gcheby(:,2:n)';
      ac(ja+(n-1)*(i-2)+(1:n-1),3*(n+1)*(m+1)+(j-1)*(n+1)+(1:n+1)) = sump;
      sumuz1 = a.*a*cheby(j,i).*ddgcheby(:,2:n)';
      sumuz2 = a*qr(:,2:n)'.*cheby(j,i).*dgcheby(:,2:n)';
      sumuz3 = 4*ddcheby(j,i).*gcheby(:,2:n)';      
      ac(ja+(n-1)*(i-2)+(1:n-1),(n+1)*(m+1)+(j-1)*(n+1)+(1:n+1)) = sumuz1+sumuz2+sumuz3;
      sumt = Ra.*cheby(j,i).*gcheby(:,2:n)';    
      ac(ja+(n-1)*(i-2)+(1:n-1),2*(n+1)*(m+1)+(j-1)*(n+1)+(1:n+1)) = sumt;      
   end
end
bc(ja+(1:(n-1)*(m-1))) = b;
ja = ja+(n-1)*(m-1);
% Tercera ecuación: divergencia
for i=2:m
   for j=1:m+1
      sumt1=a*a.*cheby(j,i).*ddgcheby(:,2:n)';
      sumt2 = a.*qr(:,2:n)'.*cheby(j,i).*dgcheby(:,2:n)';
      sumt3 = 4.*ddcheby(j,i).*gcheby(:,2:n)';
      ac(ja+(n-1)*(i-2)+(1:n-1),2*(n+1)*(m+1)+(j-1)*(n+1)+(1:n+1)) = sumt1+sumt2+sumt3;
   end
end
ja = ja+(n-1)*(m-1);
% Cuarta ecuación 
for i=2:m
   for j=1:m+1
      sumur1 = qr(:,2:n)'.*cheby(j,i).*gcheby(:,2:n)';
      sumur2 = a*cheby(j,i).*dgcheby(:,2:n)';
      ac(ja+(n-1)*(i-2)+(1:n-1),(j-1)*(n+1)+(1:n+1)) = sumur1+sumur2;
      sumuz = 2*dcheby(j,i).*gcheby(:,2:n)';    
      ac(ja+(n-1)*(i-2)+(1:n-1),(n+1)*(m+1)+(j-1)*(n+1)+(1:n+1)) = sumuz;
   end
end
ja = ja+(n-1)*(m-1); 

%    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%    %%%%%%%%%%%%%%%%Condiciones de contorno%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%      Condiciones sobre r=-1.  
%Considero Dirichlet nulas para las velocidades.
%y que condiciones sobrela presión. 
% c.c. ur=0
for i=1:m+1
   for j=1:m+1
      ac(ja+i,(j-1)*(n+1)+(1:n+1)) = cheby(j,i).*gcheby(:,1)';
   end
end
ja=ja+m+1;
% c.c uz=0
for i=1:m+1
   for j=1:m+1
      ac(ja+i,(n+1)*(m+1)+(j-1)*(n+1)+(1:n+1)) = cheby(j,i).*gcheby(1:n+1,1)';
   end
end
ja = ja+m+1;
% Bucle de elección segun diri, tomando o dirichlet o neumann. por \partial_r T = 0
% c.c. T=1
for i=1:m+1
   for J=1:m+1
      if diri==1
         ac(ja+i,2*(n+1)*(m+1)+(J-1)*(n+1)+(1:n+1))=cheby(J,i).*gcheby(1:n+1,1)';
      else
         ac(ja+i,2*(n+1)*(m+1)+(J-1)*(n+1)+(1:n+1))=cheby(J,i).*dgcheby(1:n+1,1)';
      end
   end
end
% Si no es dirichlet es cero
if diri==1
   bc(ja+(1:m+1)) = 1;
end
ja=ja+m+1;
% c.c. de Orszag sobre la presión . Considero ur identicamente cero
for i=1:m+1
   for j=1:m+1
      sump=-a*cheby(j,i).*dgcheby(:,1)';
      ac(ja+i,3*(n+1)*(m+1)+(j-1)*(n+1)+(1:n+1))=sump;
      sumur1=a.*a*cheby(j,i).*ddgcheby(:,1)';
      sumur2=qr(:,1)'.*a.*cheby(j,i).*dgcheby(:,1)';
      sumur3=4*ddcheby(j,i).*gcheby(:,1)';
      ac(ja+i,(j-1)*(n+1)+(1:n+1))=sumur1+sumur2+sumur3;
   end
end
ja=ja+m+1;
% Condiciones sobre r=1 (todos los z)
% c.c. ur=0
for i=1:m+1
   for j=1:m+1
      ac(ja+i,(j-1)*(n+1)+(1:n+1))=cheby(j,i).*gcheby(1:n+1,n+1)';
   end
end
ja=ja+m+1;

% c.c uz=0
for i=1:m+1
   for j=1:m+1
      ac(ja+i,(n+1)*(m+1)+(j-1)*(n+1)+(1:n+1))=cheby(j,i).*gcheby(1:n+1,n+1)';
   end
end
ja=ja+m+1;
% Elijo segun marque diri para la temperatura
for i=1:m+1
   for J=1:m+1
      if diri==1
         ac(ja+i,2*(n+1)*(m+1)+(J-1)*(n+1)+(1:n+1)) = cheby(J,i).*gcheby(1:n+1,n+1)';
      else
         ac(ja+i,2*(n+1)*(m+1)+(J-1)*(n+1)+(1:n+1)) = cheby(J,i).*dgcheby(1:n+1,n+1)';
      end
   end
end
ja=ja+m+1;
% condición Orszag en r=1
for i=1:m+1
   for j=1:m+1
      sump=-a*cheby(j,i).*dgcheby(:,n+1)';
      ac(ja+i,3*(n+1)*(m+1)+(j-1)*(n+1)+(1:n+1))=sump;
      sumur1=a*a*cheby(j,i).*ddgcheby(:,n+1)';
      sumur2=qr(:,n+1)'.*a*cheby(j,i).*dgcheby(:,n+1)';
      sumur3=4*ddcheby(j,i).*gcheby(:,n+1)';
      sumur4=N(n+1).*N(n+1).*cheby(j,i).*gcheby(:,n+1)';
      ac(ja+i,(j-1)*(n+1)+(1:n+1))=sumur1+sumur2+sumur3+sumur4;
   end
end
% Condición sobre la presión
ac(ja+4,:) = 0;
for j=1:m+1     
   ac(ja+4,3*(n+1)*(m+1)+(j-1)*(n+1)+(1:n+1))=cheby(j,4).*gcheby(1:n+1,n+1)';
end   
ja=ja+m+1;
% Condiciones para la capa z=-1 (todos los r menos los de los bordes
% ur=0
for j=1:m+1
   ac(ja+(1:n-1),(j-1)*(n+1)+(1:n+1))=cheby(j,1).*gcheby(:,2:n)';
end
ja=ja+n-1;
%uz=0
for j=1:m+1
   ac(ja+(1:n-1),(n+1)*(m+1)+(j-1)*(n+1)+(1:n+1))=cheby(j,1).*gcheby(:,2:n)';
end
ja=ja+n-1;
% Decrecimiento lineal en la frontera
for j=1:m+1
   ac(ja+(1:n-1),2*(n+1)*(m+1)+(j-1)*(n+1)+(1:n+1))=cheby(j,1).*gcheby(:,2:n)';
end
bc(ja+(1:n-1))=((Tmin-Tmax)*a/2.*(1./N(2:n)-a/2))/deltat+1;
ja=ja+n-1;
% Condición divergencia nula
for j=1:m+1
   sumur1=qr(:,2:n)'.*cheby(j,1).*gcheby(:,2:n)';
   sumur2=a*cheby(j,1).*dgcheby(:,2:n)';
   ac(ja+(1:n-1),(j-1)*(n+1)+(1:n+1))=sumur1+sumur2;  
   sumuz=dcheby(j,1).*gcheby(:,2:n)';    
   ac(ja+(1:n-1),(n+1)*(m+1)+(j-1)*(n+1)+(1:n+1))=sumuz;  
end
ja=ja+n-1;
% Condiciones para la capa z=1 (todos los r menos los de los bordes)

% 2*\partial_z ur+M*a*partial_r T=0
for j=1:m+1
   ac(ja+(1:n-1),(j-1)*(n+1)+(1:n+1)) = 2*dcheby(j,m+1).*gcheby(:,2:n)';
   ac(ja+(1:n-1),2*(n+1)*(m+1)+(j-1)*(n+1)+(1:n+1)) = Mar.*a.*cheby(j,m+1).*dgcheby(:,2:n)';
end
ja=ja+n-1;

% uz=0
for j=1:m+1
   ac(ja+(1:n-1),(n+1)*(m+1)+(j-1)*(n+1)+(1:n+1))=cheby(j,m+1).*gcheby(:,2:n)';
end
ja=ja+n-1;
% Condición de biot
for j=1:m+1
   sumT1 = dcheby(j,m+1).*gcheby(:,2:n)';
   sumT2 = biot*cheby(j,m+1).*gcheby(:,2:n)';
   ac(ja+(1:n-1),2*(n+1)*(m+1)+(j-1)*(n+1)+(1:n+1))=2*sumT1+sumT2;
end
ja=ja+n-1;
% Condición Orszag
for j=1:m+1
  sump=-2*dcheby(j,m+1).*gcheby(:,2:n)';
   ac(ja+(1:n-1),3*(m+1)*(n+1)+(j-1)*(n+1)+(1:n+1))=sump;
   sumuz1=a*a*cheby(j,m+1).*ddgcheby(:,2:n)';
   sumuz2=a*qr(1:n+1,2:n)'.*cheby(j,m+1).*dgcheby(:,2:n)';
   sumuz3=4*ddcheby(j,m+1).*gcheby(:,2:n)';     
   ac(ja+(1:n-1),(n+1)*(m+1)+(j-1)*(n+1)+(1:n+1))=sumuz1+sumuz2+sumuz3;
   sumt=Ra.*cheby(j,m+1).*gcheby(:,2:n)';    
   ac(ja+(1:n-1),2*(n+1)*(m+1)+(j-1)*(n+1)+(1:n+1))=sumt;        
end
bc(ja+(1:n-1))=b;
ja=ja+n-1;
%Solución del sistema y obtención de las funciones.
solucion = ac\bc;
coefur = reshape(solucion(1:(n+1)*(m+1)),n+1,m+1)';
coefuz = reshape(solucion((n+1)*(m+1)+1:2*(n+1)*(m+1)),n+1,m+1)';
coeft = reshape(solucion(2*(n+1)*(m+1)+1:3*(n+1)*(m+1)),n+1,m+1)';
coefp = reshape(solucion(3*(n+1)*(m+1)+1:4*(n+1)*(m+1)),n+1,m+1)';
% Predimensionalización de las matrices de datos
ur0 = (cheby')*coefur*gcheby;
uz0 = (cheby')*coefuz*gcheby;
t0 = (cheby')*coeft*gcheby;
p0 = (cheby')*coefp*gcheby;
%
% ---- Final function aprox2

% --- Principio function derpoly

function [A,poly]=derpoly(maximo)
% Esta función construye los polinomios de chebychev y una matriz de derivación
% Entradas
% maximo		-> Mayor grado del polinomio
% Salidas
% A			-> Matriz de derivación algebraica
% poly		-> Matriz con los coeficentes de los polinomios de chebyshev
A = zeros(maximo+1);A(1:maximo,2:maximo+1) = diag(maximo:-1:1);
poly = zeros(maximo+1);C = zeros(maximo+2);
poly(1,:) = [zeros(1,maximo),1];
poly(2,:) = [zeros(1,maximo-1),1,0];
for i = 1:maximo-1;
   C(i+2,:) = 2.*[poly(i+1,:),0]-[0,poly(i,:)];
   poly(i+2,:) = [C(i+2,2:maximo+2)];
end

% --- Final function derpoly 
%
%% ---- Principio function derevol. Función técnica
function [varargout] = derevol(a,C)
% Esta función sustituye a deralgornm. Obtiene las derivadas de una función cuyta
% malla está en los puntos de chebyshev
% Entradas %
% a		-> Matriz de la que se quieren obtener las derivadas
% gc 		-> Matriz de colocación
% cheby 	-> Tensores con los puntos de colocación.
% Es un algoritmo un poco rarito
% Parte comun del algoritmo
r = C{1};
z = C{2} ;
cheby = C{3};
dcheby = C{4};
ddcheby = C{5};
gcheby = C{6};
dgcheby = C{7};
ddgcheby = C{8};
icheby = C{9};
igcheby = C{10};
[M,N] = size(a);
x = icheby*a*igcheby;
n = N-1;
m = M-1;
if nargout == 2				% Esto significa que es la presión
   varargout(1) = {(cheby')*x*dgcheby};
   varargout(2) = {(dcheby')*x*gcheby};
else
	auxr = (cheby')*x*dgcheby;
	auxz = (dcheby')*x*gcheby;
	varargout(1) = {auxr};  			% \partial_r
	varargout(2) = {auxz};				% \partial_z
	varargout(3) = {(cheby')*x*ddgcheby};% \partial_rr
	varargout(4) = {(ddcheby')*x*gcheby};%\partial_zz
	Tra = zeros(N,N,M);
	Tza = zeros(N,N,M);
	for i=1:M
		Tra(:,:,i) = auxr(ones(N,1),:);
		Tza(:,:,i) = auxz(ones(N,1),:);
	end
	varargout(5) = {Tra};
	varargout(6) = {Tza};
end

% --- Final function derevol ---
%
% --- datos fijos del problema ----
function C=datos(n,m)
maximo = max(n,m);
% Puntos de colocación de Gauss-Lobato, y funciones técnicas. 
r = [-1 cos(pi.*(n-1:-1:1)./n) 1];
z = [-1 cos(pi.*(m-1:-1:1)./m) 1];
hi = 1;
% polinomios de tchebicheff y matriz de derivación
[A,poly] = derpoly(maximo);
% Predimensionalización de las matrices de datos
dpoly = poly*A;
ddpoly = dpoly*A;
cheby = zeros(m+1);
dcheby = zeros(m+1);
ddcheby = zeros(m+1);
gcheby = zeros(n+1);
dgcheby = zeros(n+1);
ddgcheby = zeros(n+1);
% Construimos las evaluaciones en los puntos de colocación.
for i = 1:m+1
   cheby(i,:) = [polyval(poly(i,:),z)];
   dcheby(i,:) = [polyval(dpoly(i,:),z)];
   ddcheby(i,:) = [polyval(ddpoly(i,:),z)];
end
for i = 1:n+1
   gcheby(i,:) = [polyval(poly(i,:),r)];
   dgcheby(i,:) = [polyval(dpoly(i,:),r)];
   ddgcheby(i,:) = [polyval(ddpoly(i,:),r)];
end
% Construci´on del Cell para almacenar
icheby = inv(cheby');
igcheby = inv(gcheby);
C = cell(1,10);
C{1} = r;
C{2} = z;
C{3}= cheby;
C{4}= dcheby;
C{5}= ddcheby;
C{6}= gcheby;
C{7}= dgcheby;
C{8}= ddgcheby;
C{9}= icheby;
C{10} = igcheby;