[approx,exac,err]=odegalerkin(poly,cc,n)

% ODEGALERKIN Ordinary Differential Equation Solver through Galerkin
% Method.
% [APPROX,EXAC,ERR] = ODEGALERKIN(POLY,CC,N) solves Ordinary Differential 
% Equations (ODE)through Galerkin method, by inserting the characteristic 
% polynomial matrix "POLY", boundary conditions "CC" and the finite 
% quantity of approximative base functions "N".
%
% Outputs of the program are the approximative solution "APPROX", the
% analitic solution "EXAC" and the percentage error "ERR" (%).
%
% NOTE: IF BOUNDARY CONDITIONS ARE EXPRESSED IN TERMS OF DERIVATIVES, THE
% "CC" MATRIX MUST BE OF SIZE 4X4 INSTEAD OF 1X4, SIMPLY BY PUTTING ";"
% BETWEEN THE SECOND AND THIRD ELEMENT.
%
% 
% EXAMPLE 1:
% Find an approximative solution to the following ODE. Use 17 base
% functions.
% 
% 1*D2y(x)+2*Dy(x)-3*y(x)=6      y(2)=4     y(0)=-7
%
% ======Program======
% [approx]=odegalerkin([1 2 -3 6],[2 4 0 -7],17)
%
% EXAMPLE 2:
% Find an approximative solution to the following ODE. Use 25 base
% functions. Also calculate the exact solution.
%
% 1*D2y(x)+0*Dy(x)+3*y(x)=x^2+2       y(2)=4     y(0)=7
%
% ======Program======
% syms x
% [approx,exac]=odegalerkin([1 0 3 x^2+2],[2 4 0 7],25)
%
% EXAMPLE 3:
% Find an approximative solution to the following ODE. Use 4 base
% functions. Also calculate the exact solution and the percentage error.
%
% -x*D2y(x)+4*Dy(x)+3*x*y(x)=8*sin(x)      y(2)=0     y'(5)=-2
%
% ======Program======
% syms x
% [approx,exac,error]=odegalerkin([-x 4 3*x 8*sin(x)],[2 0; 5 -2],4)
%
% (c)2008, Marcos Cesar Ruggeri
% [e-mail]
% Universidad Tecnologica Nacional
% Facultad Regional Haedo
% Departamento Aeronautica

function [approx,exac,err]=odegalerkin(poly,cc,n)

syms x;
alpha=poly(1);
beta=poly(2);
gamma=poly(3);
g=poly(4);
if length(cc)==4
    x1=cc(1);
    phi1=cc(2);
    x2=cc(3);
    phi2=cc(4);
    phi=(phi2-phi1)/(x2-x1)*(x-x1)+phi1;   
elseif length(cc)==2
    cc=cc';
    x1=cc(1);
    phi1=cc(2);
    x2=cc(3);
    phiprima2=cc(4);
    phi=phiprima2*x+(phi1-phiprima2*x1);
    phi2=phiprima2*x2+(phi1-phiprima2*x1);
else error('The size of the boundary conditions matrix is not valid.')
end

phider1=diff(phi);
phider2=diff(diff(phi));
if lt(x1,x2)
    x_inf=x1;
    x_sup=x2;
    y_inf=phi1;
    y_sup=phi2;
else
    x_inf=x2;
    x_sup=x1;
    y_inf=phi2;
    y_sup=phi1;
end
for i=1:n, 
    psi(1,i)=x^(i-1); 
end
if length(cc)==4
    psi=(x-x1)*(x-x2)*psi;
elseif length(cc)==2
    psi=(x-x1)*((x-x2)^2).*psi;
else error('The size of the boundary conditions matrix is not valid.')
end
psider1=diff(psi);
psider2=diff(diff(psi));
psitrans=transpose(psi);
E1=alpha*psider2+beta*psider1+gamma*psi;
E2=alpha*phider2+beta*phider1+gamma*phi-g;
c=inv(int(psitrans*E1,x,x1,x2))*-int(psitrans*E2,x,x1,x2);
e=E1*c+E2;

fun_approx=phi+psi*c;
approx=fun_approx;
fun_approx2=inline(fun_approx);
x=(x_inf-10):.2:(x_sup+10);
y_approx=subs(fun_approx);

clear x;
syms x;
if length(cc)==4
    fun_exac=dsolve('alpha*D2y+beta*Dy+gamma*y=g','y(x1)=phi1','y(x2)=phi2','x');
elseif length(cc)==2
    fun_exac=dsolve('alpha*D2y+beta*Dy+gamma*y=g','y(x1)=phi1','Dy(x2)=phiprima2','x');
else error('The size of the boundary conditions matrix is not valid.')
end
fun_exac=subs(fun_exac);
exac=fun_exac;
fun_exac2=inline(fun_exac);
x=(x_inf-10):.2:(x_sup+10);
y_exac=subs(fun_exac);

e_norm2=sqrt(int(e^2,x1,x2));
f_norm2=sqrt(int(fun_approx^2,x1,x2));
err=abs(vpa(e_norm2/f_norm2*100,2));

plot(x,y_approx,'-o','Color','r'),grid
axis([x_inf-pi x_sup+pi -100 100])
hold on
title('\fontsize{14} Solutions of the Ordinary Differential Equation through Galerkin Method')
xlabel('x')
ylabel('\Phi(x)')
plot(x,y_exac),grid
legend('Approximative solution','Exact solution')
hold off
grid on