% OCW Universidad de Cantabria
% Cálculo Simbólico y Numérico en Ecuaciones Diferenciales
% © Cursos 2001-2014,  M. Eugenia Pérez, UC
%
function y = elementosfinitos(n);
%
h=1/(n+1);
%
x=ones(n+3,1);
for i=1:n+2
   x(i)=(i-1)*h;
end
nodos=x(1:n+2)
%
syms t; 
q=-1; 
r=t^2; 
R=ones(n+1,1); 
S=ones(n+1,1);
Q=ones(n+1,1);
f=ones(n+1,1);
a1=ones(n+1,1);
amenos1=ones((n),1);
mas1=ones((n),1);
%
for i=1:n+1
   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
f(n+1)=double((1/h)*int(r*(t-x(n+1)),x(n+1),x(i+2)));
S1menos=S(2:n+1);
%
a1= (1/h^2)*(-2*h*ones(n+1,1)+Q+R);  
amenos1=(1/h^2)*(h*ones((n),1)-S1menos);  
amas1=amenos1; 
%
%
A=diag(a1)+diag(amas1,1)+diag(amenos1,-1);  
A(n+1,n+1)=(1/h^2)*(-h+Q(n+1));
%
matriz_sistema=A
termino_independiente=f
%
c=A\f;
coeficientes_solucion=c
v=dsolve('D2y-y=t^2','y(0)=0','Dy(1)=0');
% 
solu_exacta=v
% 
%
ezplot(v,[0,1])
%
hold on
plot(nodos,[0;c],'r.')
hold off
%
solucionaproximadaennodos=[0;c]
solucionexactaennodos=subs(v,t,nodos)
%
error=solucionaproximadaennodos-solucionexactaennodos
 
