function y = elementosfinitos(n)
%CALCULO SIMBOLICO Y NUMERICO EN E.D.
%%V.MATLAB 6.5
%y = felementosfinitos(n)
%
%calculo de la aproximacion de la solucion del problema de contorno 
%%y''-y=x^2, x\in(0,1); y(0)=0,y(1)=0
%utilizando el método de los elementos finitos: se toman las funciones
%de base phi_i(x) que son continuas en [0,1], se anulan en x=0 y x=1
%son polinomios de grado 1 en los subintervalos [x_i,x_{i+1}], con x_i los nodos 
%y toman el valor 1 en el nodo x_i y 0 en el resto ''hat functions''
%
%lee n el número de nodos del intervalo (0,1)
%
%calcula la matriz del sistema A y el termino independiente f y 
%nos devuelve un vector c que es la solucion aproximada de y(x) en los nodos x_i=i*h
% 
%definiendo los nodos
h=1/(n+1);
x=ones(n+2,1);
for i=1:n+2
   x(i)=(i-1)*h;
end
%
%Calculando la matriz de coeficientes A y el término independiente:
%
syms t; %variable simbolica
q=-1;%termino que aparece multiplicando a y(x): y''+q(x)y=r(x)
r=t^2;%termino independiente: y''+q(x)y=r(x)
R=ones(n,1);%inicializando vectores
S=ones(n,1);
Q=ones(n,1);
f=ones(n,1);
a1=ones(n,1);
amenos1=ones((n-1),1);
mas1=ones((n-1),1);
%definiendo las integrales sobre cada segmento [x_i,x_{i+1}]
for i=1:n
   R(i)=double(int(q*(t-x(i+2))^2,x(i+1),x(i+2)));
   Q(i)= double(int(q*(t-x(i))^2,x(i),x(i+1)));
   S(i)=double(int(q*(t-x(i))*(t-x(i+1)),x(i),x(i+1)));
   f(i)=double((1/h)*int(r*(t-x(i)),x(i),x(i+1))+(1/h)*int(r*(-t+x(i+2)),x(i+1),x(i+2)));
end
%S1mas=S(1:n-1);
S1menos=S(2:n);
%
%definiendo la matriz tridiagonal A:  vectores que definen  las diagonales
%diagonal principal
a1= (1/h^2)*(-2*h*ones(n,1)+Q+R); %A(i,i)
%subdiagonal inferior
amenos1=(1/h^2)*(h*ones((n-1),1)-S1menos); %A(i,i-1)
%subdiagonal superior
amas1=amenos1;  %A(i,i+1)
%
%
A=diag(a1)+diag(amas1,1)+diag(amenos1,-1); %matriz tridiagonal A
%
matriz_sistema=A
termino_independiente=f
%
%resolviendo el sistema
%
c=A\f;
coeficientes_solucion=c
%
Nodo_y_Aproximacionsolucion=[x,[0;c;0]]
%
%definiendo la solucion exacta
v=dsolve('D2y-y=t^2','y(0)=0','y(1)=0');
solu_exacta=v  %solucion exacta v
%dibujando las dos soluciones
%
ezplot(v,[0,1])
%
hold on
%
plot(x,[0;c;0],'--')
hold off
%Comparando soluciones
%
nodos=x
solucionaproximadaennodos=[0;c;0];
solucionexactaennodos=subs(v,t,x);
%
%error cometido
%
error=solucionaproximadaennodos-solucionexactaennodos
%double(error)
%cambio de formato en la salida de datos
%%format long 
%%nodo_error=[nodos,error]
%%OBSERVACIONES:
%1- La matriz A es simetrica; A es tridiagonal. El calculo de las integrales en [0,1] 
%   se lleva al calculo de integrales sobre cada segmento [x_i,x_{i+1}].
%   Vale con calcular la diagonal por encima de la diagonal principal
%
%2- En este ejemplo concreto podemos calcular explicitamente las integrales.
%  Tambien, en este ejemplo concreto, mediante resolucion, se pueden comparar 
%  la solucion exacta y(x) con la aproximada u=sum_{i=1}^n c_i\phi_i(x) 
%
%3- El vector solucion c nos da la solucion aproximada de y(x) en los nodos x_i=i*h
%%%%mas concretamente, c(i) nos da la solucion aproximada en x(i+1); i=1,2,..,n
%%%%c(1)=y(0)=0; c(n+2)=y(1)=0

%
%4-Con muy pocas modificaciones se podria resolver el problema mas general
%%(p(x)y')'+q(x)y=f(x), x\in(a,b), y(a)=0,y(b)=0; 
%%Asi como otras condiciones de contorno: y'(a)=0, o y(b)-y'(b)=0,....

