function [uh] = solveP_FEM(caso, hmax, dx, dy)
	O = ... %definimos un dominio inicial
	A = 1;
	rhs=0;

	switch (caso)
		case 1
			[P,T,ND,B,O] = ... % obtener triangulacion P, T, ND, B (puntos, triangulos, nodos dirichlet, bases) tal que h < hmax. Ademas, definir el dominio omega (O) corresopondiente a partir de los puntos que configuran el borde, si fuera necesario.
			[A,rhs] = fem(P, T, ND, B, @f1);

		case 2
			[P,T,ND,B,O] = ... % obtener triangulacion P, T, ND, B (puntos, triangulos, nodos dirichlet, bases) tal que h < hmax. Ademas, definir el dominio omega (O) corresopondiente a partir de los puntos que configuran el borde, si fuera necesario.
			[A,rhs] = fem(P, T, ND, B, @f2);

		case 3
			[P,T,ND,B, O] = ... % obtener triangulacion P, T, ND, B (puntos, triangulos, nodos dirichlet, bases) tal que h < hmax. Ademas, definir el dominio omega (O) corresopondiente a partir de los puntos que configuran el borde, si fuera necesario.
			[A,rhs] = fem(P, T, ND, B, @f3);

		otherwise
			uh=0;
	end

	uh=A\rhs
	... % graficar la solucion
end

function [f]=f1(x,y)
	... %definicion de la densidad de carga correspondiente al 
		%caso 1
end

function [f]=f2(x,y)
	... %definicion de la densidad de carga correspondiente al 
		%caso 2
end

function [f]=f3(x,y)
	... %definicion de la densidad de carga correspondiente al 
		%caso 3
end
