Daubechies Polynomial, Scaling Function and Wavelet

First define the polynomial
> p0:=(1+sqrt(3))/4; p1:=(3+sqrt(3))/4; p2:=(3-sqrt(3))/4; p3:=(1-sqrt(3))/4;

                       p0 := 1/4 + 1/4 sqrt(3)


                       p1 := 3/4 + 1/4 sqrt(3)


                       p2 := 3/4 - 1/4 sqrt(3)


                       p3 := 1/4 - 1/4 sqrt(3)

> P:=z->(p0+p1*z+p2*z^2+p3*z^3)/2;

                                               2           3
         P := z -> 1/2 p0 + 1/2 p1 z + 1/2 p2 z  + 1/2 p3 z

Now show that P(z) satisfies the key conditions
needed so that the scaling function and its
translates are orthonormal.
> z:=exp(I*t); zc:=exp(-I*t);

                            z := exp(I t)


                           zc := exp(-I t)

> expr:=P(z)*P(zc) +P(-z)*P(-zc):
> simplify(%);

                                  1

> f:=(Re(P(z)))^2+(Im(P(z)))^2;

  f := (1/8 + 1/8 sqrt(3) + 1/2 Re((3/4 + 1/4 sqrt(3)) exp(I t)

                                       2
         + (3/4 - 1/4 sqrt(3)) exp(I t)

                                       3  2
         + (1/4 - 1/4 sqrt(3)) exp(I t) ))  + 1/4 Im(

                                                                   2
        (3/4 + 1/4 sqrt(3)) exp(I t) + (3/4 - 1/4 sqrt(3)) exp(I t)

                                       3 2
         + (1/4 - 1/4 sqrt(3)) exp(I t) )

> plot(f,t=-Pi/2..Pi/2,y=0..5);

Here is the iterative process for constructing
the scaling function. Start with phi_0 = Haar scaling function
> phi_0:=x->piecewise(x<0,0,x<1,1,0);

            phi_0 := x -> piecewise(x < 0, 0, x < 1, 1, 0)

> plot(phi_0,-1..3);

Now iterate the scaling equation to define
phi_1, phi_2 etc.
> phi_1:=x->p0*phi_0(2*x)+p1*phi_0(2*x-1)+p2*phi_0(2*x-2)+p3*phi_0(2*x-3);

  phi_1 := x -> p0 phi_0(2 x) + p1 phi_0(2 x - 1) + p2 phi_0(2 x - 2)

         + p3 phi_0(2 x - 3)

> plot(phi_1,-1..3, numpoints=100);

Here is the wavelet function.
> psi:=x->-p0*phi_5(2*x-1)+p1*phi_5(2*x)-p2*phi_5(2*x+1)+p3*phi_5(2*x+2);
>   psi := x -> -p0 phi_1(2 x - 1) + p1 phi_1(2 x) - p2 phi_1(2 x + 1)
> 
>          + p3 phi_1(2 x + 2)
> 

  psi := x -> -p0 phi_5(2 x - 1) + p1 phi_5(2 x) - p2 phi_5(2 x + 1)

         + p3 phi_5(2 x + 2)

> plot(psi,-5..8, numpoints=100, tickmarks=[0,0]);
>