function y= galerkin(n)
%%CALCULO SIMBOLICO Y NUMERICO EN E.D.
%V.MATLAB 6.5
%function c = galerkin(n)
%calcula la aproximacion del problema de contorno
%y''- y=x^2, x\in(0,1), y(0)=0,y(1)=0. 
%por el metodo de Galerkin, 
%utilizando como funciones de base las funciones propias del sistema
%%y''+\lambda y=0, x\in(0,1), y(0)=0,y(1)=0. 
%
%Estas explicitamente son   sin(k*pi*x), k=1,2,..
%(asociadas a los valores propios lambda=k^2pi^2)
%
%La funcion lee n el numero de
%elementos de la base y nos devuelva el vector c que nos permite
%calcular la aproximacion numerica de la solucion;
%
%calcula la solucion aproximada y  la compara con la exacta.
%
%Hay que resolver el sistema: Ac=f, 
%donde A es la matriz de coeficientes; f el dato del segundo miembro
%la aproximacion de la solucion  se obtiene: u=y_n(t)=sum_{j=1}^n (c_j*sin(j*pi*t))
%
%%definiendo la matriz A y el vector columna f
%
syms t %variable simbolica
A=ones(n,n); %incicializando vectores y matrices
f=ones(n,1);
for i=1:n
   f(i)=-double(int(sin(i*pi*t)*t^2,0,1));
   for j=1:n
      if i==j
         A(i,j)=double(int((i*pi*cos(i*pi*t))^2,0,1)+int((sin(pi*i*t))^2,0,1));
               else
         A(i,j)=double(int(i*pi*cos(i*pi*t)*j*pi*cos(j*pi*t),0,1)+ int(sin(i*pi*t)*sin(j*pi*t),0,1));
       end
      end
   end
 %
matriz_sistema=A
termino_independiente=f
%
%resolviendo el sistema
%
c=A\f;
coeficientes_solucion=c
%
% definiendo la solucion aproximada
%
u=0;
for i=1:n
   u=u+c(i)*sin(pi*i*t);
end
%
solu_aproximada=vpa(u,4) %u solucion aproximada (4 digitos solo)
%
%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
%
dt=0:0.05:1;
du=subs(u,t,dt);
ddu=double(du);
plot(dt,ddu,'--')%%plot(dt,du,'--') directamente no funciona en la version 5.0
%
hold on  
%ezplot(u,[0,1])%tambien dibuja la solucion u
ezplot(v,[0,1])
%
hold off
%
%%comparando los coeficientes c(i) con los de Fourier alpha(i) de la solucion exacta v
%
alpha=ones(n,1);
for i=1:n
   alpha(i)=double(int(sin(i*pi*t)*(-2-t^2-(-3+2*cosh(1))/sinh(1)*sinh(t)+2*cosh(t)),0,1)/int((sin(pi*i*t))^2,0,1));
end
coeficientes_fourier=alpha
%
% error cometido al tomar los coeficientes c(i) por alpha(i)
% 
errorfourier_galerkin=alpha-c
%
%%OBSERVACIONES:
%1- La matriz A es simetrica; por ser los coeficientes constantes y tener una 
%   base ortogonal en (0,1), en este caso, la matriz del sistema es diagonal 
%
%2- En este ejemplo concreto podemos comparara con la solucion exacta
%
%3-Con muy pocas modificaciones se podria resolver el problema más general
%%(p(x)y')'+q(x)y=f(x), x\in(a,b), y(a)=0,y(b)=0. 
%
