kill(all)$
load(functs)$

prepare:truebasic code (implementation),to do integration with simple rule
        for quadratic element in one dimension.
rule : a*x^n -> integration for x = ' a*x^(n+1)/(n+1) '
(need to find coefficient of a,see 'powerseries above')
peter.vlasschaert@gmail.com ,16/07/2016

y1:a1*y^2+a2*y+a3$   /* mu(y) */
                     /* see : quadratic element (website) */
y2:a4*y+a5$          /* quadratic element (element matrix),derivatives 'N1_D,N2_D,N,N3_D'*/
y3:a6*y+a7$          /* quadratic element (element matrix),derivatives 'N1_D,N2_D,N,N3_D'*/;

yy:y1*y2*y3;

yyy:ratsimp(yy,y);

b0:ev(yyy,y=0);
b1:coeff(yyy,y);
b2:coeff(yyy,y^2);
b3:coeff(yyy,y^3);
b4:coeff(yyy,y^4);

yyyy:integrate(yyy,y);

c1:coeff(yyyy,y^1);
c2:coeff(yyyy,y^2);
c3:coeff(yyyy,y^3);
c4:coeff(yyyy,y^4);
c5:coeff(yyyy,y^5);

integration between two values in the interval : 'beginy = value_1 < y < endy =value_2'

yyyyy:integrate(yyy,y,beginy,endy);